Report

周波数領域における画像処理 林沼 勝利 廣安 知之 山本 詩子 2014 年 6 月 4 日 IS Report No. 2014060401 Report
Medical Information System Labratry Abstract
画像処理において,指紋認証や虹彩認識などの生体認証,テクスチャのパターン認識などは画像の
周波数領域における情報が重要となる.本稿では画像の空間周波数を解析する手法としてのフーリエ
変換やウェーブレット変換について述べる.
目次
第 1 章 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
第 2 章 フーリエ変換
. . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.1
離散フーリエ変換 . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.2
実装例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.3
周波数フィルタリング . . . . . . . . . . . . . . . . . . . . . . .
5
2.3.1
ローパスフィルタ . . . . . . . . . . . . . . . . . . . . . . . .
5
2.3.2
ハイパスフィルタ . . . . . . . . . . . . . . . . . . . . . . . .
6
第 3 章 ウェーブレット変換 . . . . . . . . . . . . . . . . . . . . . . . .
7
3.1
短時間フーリエ変換 . . . . . . . . . . . . . . . . . . . . . . . .
7
3.2
連続ウェーブレット変換 . . . . . . . . . . . . . . . . . . . . . .
7
3.3
画像処理におけるウェーブレット変換 . . . . . . . . . . . . . . . . .
7
3.3.1
ハールウェーブレット . . . . . . . . . . . . . . . . . . . . . .
7
3.3.2
ガボールウェーブレット . . . . . . . . . . . . . . . . . . . . .
9
第 1 章 はじめに
音や電波などの信号の周波数とは,一定時間内に変化する振動の回数のことを指す.画像において
も同様の考え方が可能であり,画像の濃度値をパラメータとしてグラフを書くことにより,明るさの
変化に応じて波を形成することができる (Fig. 1.0.1).そこで,
「一定の長さで,どのくらい明るさの
強弱が変化するのか」で画像の周波数を決定することができる.画像では,時間を空間に置き換えて
周波数を定義する.これを空間周波数という.本稿では,この空間周波数の解析を行う手法としての
フーリエ変換とウェーブレット変換について述べる.
なお,本稿の第 2 章は「ディジタル画像処理」(CG-ARTS 協会)1) および「C 言語で学ぶ医用画像
処理」(オーム社)2) ,第 3 章は「やり直しのための通信数学」(CQ 出版)3) を参考に作成している.
(a) 2 次元濃度パターン
(b) 横一列の濃度値変化
Fig. 1.0.1 sin 波の画像の変化 (自作)
2
第 2 章 フーリエ変換
2.1
離散フーリエ変換
画像が持っている空間周波数を利用した処理をするためには,画像をまず空間周波数領域へ変換し,
処理を行った後に変換して,元の空間領域の画像に戻す.この空間領域と空間周波数領域の相互変換
に用いるのがフーリエ変換とフーリエ逆変換のフーリエ変換対である.まず,1 次元のフーリエ変換
対は次の式であらわされる.
∫
∞
F (u) =
f (x)e−j2πux dx
(2.1)
F (u)ej2πux du
(2.2)
−∞
∫ ∞
f (x) =
−∞
これらを 2 次元に拡張した 2 次元フーリエ変換対は以下の式であらわされる.
∫ ∞∫ ∞
f (x, y)e−j2π(ux+vy) dxdy
F (u, v) =
−∞
∫
∞
(2.3)
−∞
∫
∞
F (u, v)ej2π(ux+vy) dudv
f (x, y) =
(2.4)
−∞
−∞
これらのフーリエ変換は連続信号を対象としているため,離散信号を扱う場合には場合には離散フー
リエ変換(Discrete Fourier Transform: DFT)が用いられる.1 次元の DFT は以下の式であらわさ
れる.
F (u) =
N
−1
∑
f (x)e−j
2πux
N
(2.5)
x=0
f (x) =
N −1
2πux
1 ∑
F (u)ej M
N
(2.6)
u=0
また,2 次元の DFT は以下の式であらわされる.
F (u, v) =
N
−1 M
−1
∑
∑
f (x, y)e−j
2πux
M
e−j
2πvy
N
(2.7)
y=0 x=0
f (x, y) =
N −1 M −1
2πvy
2πux
1 ∑ ∑
F (u, v)ej M ej N
MN
(2.8)
v=0 u=0
ここで,画像に対して DFT を行う場合を考える.式 (2.7) について
fˆ(u, v) =
M
−1
∑
f (x, y)e−j
2πux
M
(2.9)
2πvy
fˆ(u, v)e−j N
(2.10)
x=0
F (u, v) =
N
−1
∑
y=0
3
2.2 実装例
第 2 章 フーリエ変換
と表すことができる.これは,画像に対してまず x 方向について 1 次元 DFT を行い,その結果に対
して y 方向に 1 次元 DFT を行えばよいことがわかる.式 (2.9) についてオイラーの公式を用いるこ
とにより
fˆ(u, v) =
M
−1
∑
f (x, y)(cos(
x=0
2πux
2πvy
) − j sin(
))
M
N
(2.11)
と展開することができる.y 方向についても同様である.これにより積和演算として DFT を計算す
ることができる.これらを用いてフーリエ変換の結果を画像として表示する場合,次の式であらわさ
れるパワースペクトルを用いて表示する.
|F (u, v)|2 = Re {F (u, v)}2 + Im {F (u, v)}2
(2.12)
通常,このパワースペクトルを画像として表示する場合,対数にして表示する.また,2D-DFT の結
果は四隅に DC 成分が来るため,Fig. 2.1 に示すように中央に DC 成分がくるようにデータを入れ替
える.
低
高高
低
低
高
高
低
Fig. 2.1.1 周波数領域の入れ替え(自作)
2.2
実装例
Python を用いて画像にフーリエ変換をかけるコードを以下に示す.なお,DFT は O(n2 ) の計算量
が必要となるため,ここでは DFT を高速化した高速フーリエ変換(Fast Fourier Transform: FFT)
を用いる.
ソースコード 2.1 FFT
1
2
3
# -* - coding : utf -8 -* from PIL import Image
import numpy as np
4
5
6
7
8
9
10
11
12
13
14
im = Image . open ( " lena . jpg " )
im = np . array ( im )
ft = np . fft . fft2 ( im )
# 2 次 元 FFT
ft = np . fft . fftshift ( ft )
#DC成分が中央にくるようにデータを入れ替え
Pow = np . abs ( ft )**2
# パワースペクトルを計算
Pow = np . log10 ( Pow )
# 対数表示
Pmax = np . max ( Pow )
Pow = Pow / Pmax *255
pow_img = Image . fromarray ( np . uint8 ( Pow ))
pow_img . save ( ’ lena_fft . jpg ’)
この実行結果を Fig. 2.2.1 に示す.
4
2.3 周波数フィルタリング
第 2 章 フーリエ変換
(a) 入力画像 (lena.jpg)
(b) 出力画像 (lena fft.jpg)
Fig. 2.2.1 ソースコード 2.1 の結果(自作)
2.3
周波数フィルタリング
フーリエ変換の各周波数成分は,画像に含まれるそれぞれの空間周波数の強さを表している.そこ
で,フーリエ変換後の各周波数成分の大きさを各成分ごとに変えることにより,元の画像の性質を変
化させることができる.このような処理を周波数フィルタリングという.元の画像をフーリエ変換し
たものを F (u, v),フィルタリングの出力を G(u, v) とするとき,周波数フィルタリングは以下の式で
表される.
G(u, v) = F (u, v) · H(u, v)
(2.13)
ここで,H(u, v) が周波数フィルタで,各周波数成分にかけられる倍率に相当する.
2.3.1
ローパスフィルタ
画像に含まれる空間周波数成分のうち,高周波成分を遮断して低周波成分のみを残すようなフィル
タをローパスフィルタ(lowpass filter: LPF)という.Fig. 2.3.1 にその一例を示す.Fig. 2.3.1 (a)
はフィルタの関数を 3 次元的にプロットしたものである.Fig. 2.3.1 (b) はレナの画像に対して FFT
をかけ,LPF を通し,IFFT を行った結果である.高周波成分が取り除かれ,ゆるやかな変動だけが
残り,平滑化された画像が得られていることが分かる.
(a) LPF フィルタ関数のプロット
(b) LPF の結果
Fig. 2.3.1 LPF の例(自作)
5
2.3 周波数フィルタリング
2.3.2
第 2 章 フーリエ変換
ハイパスフィルタ
ローパスフィルタとは逆に,画像に含まれる空間周波数成分のうち,低周波成分を遮断して高周波
成分のみを残すようなフィルタをハイパスフィルタ(highpass filter: HPF)という.Fig. 2.3.2 にそ
の一例を示す.Fig. 2.3.2 (a) はフィルタ関数を 3 次元的にプロットしたものである.Fig. 2.3.2 (b)
はレナの画像に対して FFT をかけ,HPF を通し,IFFT を行った結果である.低周波成分が取り除
かれ,細かな濃度変化だけが残り,エッジが強調された画像が得られていることが分かる.
(a) HPF フィルタ関数のプロット
(b) HPF の結果
Fig. 2.3.2 HPF の例(自作)
6
第 3 章 ウェーブレット変換
3.1
短時間フーリエ変換
1 次元の信号に対してフーリエ変換を行った場合,周波数空間では時間情報が失われている.そこ
で,信号の一部を窓関数を用いて切りだし,この窓をずらすことによりスペクトルの時間変化を求め
ることができる.これを短時間フーリエ変換(Short-Time Fourier Transform: SIFT)という.窓関
数を ω(t) とすると SIFT は次の式で表現できる.
∫ ∞
F (b, ω) =
ω(t − b)f (t)e−jωt dt
(3.1)
−∞
しかし,SIFT は解像度が限られるという問題点がある.幅の広い窓は周波数分解能が高いが,時間
分解能は低くなってしまう.逆に幅の狭い窓は時間分解能は高くなるが,周波数分解脳が低くなって
しまう.これを不確定性原理という.
これは画像のフーリエ変換においても同様のことが言える.画像のフーリエ変換の結果は位置の
情報が失われてしまっている.そこで,画像を窓関数により局所化することにより位置情報を残して
フーリエ変換をすることができる.しかし,1 次元の信号の SIFT と同様に小さい窓を使用すると周
波数が求まらなくなり,大きな窓を使用すると場所を定めることができなくなる.
3.2
連続ウェーブレット変換
前節の問題点を解決するため,ウェーブレット変換が考えられるようになった.平均値が0で,原
点 t = 0 の周りに局在する関数を ψ(t) として,信号データ x(t) との内積を
∫ ∞
1 ¯ x−b
√ ψ(
(Wψ f )(b, a) =
)f (t)dt
a
|a|
−∞
(3.2)
で定義し,これを連続ウェーブレット変換と呼ぶ.ここで,ψ¯ は ψ の複素共役を表し,ψ(t) をアナ
ライジング・ウェーブレットまたは単にウェーブレットと呼ぶ.パラメータ b はウェーブレットの時
間シフトを表し,a はウェーブレットを拡大・縮小する倍率を表す.つまり,ウェーブレット変換は
ψ(t) を拡大・縮小した関数を用いて信号の性質を表している.
3.3
3.3.1
画像処理におけるウェーブレット変換
ハールウェーブレット
アナライジングウェーブレットに Haar 関数 (式 (3.3)) を用いたウェーブレットをハールウェーブ
レット(Haar Wavelet)と呼ぶ.



1
0 ≦ t < 1/2



ψ(x) = −1 1/2 ≦ t < 1




0
otherwise
7
(3.3)
3.3 画像処理におけるウェーブレット変換
第 3 章 ウェーブレット変換
Haar 関数をグラフで表すと Fig. 3.3.1 となる.Haar ウェーブレット変換は平均と平均差分を求める
計算となる.このとき,平均は低周波数成分を,平均差分は高周波成分を取り出すフィルタの働きに
等価であることが知られている.画像に対して行う場合,フーリエ変換と同様 x 方向に 1 次元ウェー
ブレット変換を行った後,その結果に対して y 方向にウェーブレット変換を行う.この様子を Fig.
3.3.2 に示す.一度ウェーブレット変換を行うと HL には縦方向のエッジ,LH には横方向のエッジ,
HH には対角方向のエッジが得られる.LL 成分は原画像を 1/2 に縮小した画像となる.また,LL 成
分に対してウェーブレット変換を行うと,さらに 4 分割された結果が得られる.これにより,1 枚の
画像を複数の解像度を有する画像成分の重ね合わせにより,階層的に表現することができる.このよ
うに,解像度を変えて信号を解析することを多重解像度解析という.Fig. 3.3.3 にレナの画像に対し
て 2 階層 Haar ウェーブレット変換を行った結果を示す.
Fig. 3.3.1 Haar 関数(自作)
Fig. 3.3.2 画像に対する Haar ウェーブレット変換(自作)
Fig. 3.3.3 Haar ウェーブレットの実行結果(自作)
8
3.3 画像処理におけるウェーブレット変換
3.3.2
第 3 章 ウェーブレット変換
ガボールウェーブレット
アナライジングウェーブレットに Gabor 関数 (式 (3.4)) を用いたウェーブレットをガボールウェー
ブレット(Gabor Wavelet)と呼ぶ.
t2
1
ψ(x) = √ 2 e− 2σ2 e−jωx
2 πσ
(3.4)
式 (3.4) で σ = 8 とした時の実部のグラフを Fig. 3.3.4 に示す.Gabor 関数はガウス関数と正弦波の
積となっている.画像処理においてガボールウェーブレットを行う場合,以下に示す式をよりフィル
タバンクを作成する.
)
(
)
( ′2
x′
x + γ 2 y ′2
g(x, y; λ, θ, ψ, σ, γ) = exp −
cos 2π + ψ
2σ 2
λ
ここで
 
  
′
cos θ − sin θ
x
x
 
 =
y′
sin θ cos θ
y
(3.5)
(3.6)
である.λ は波長のコサイン成分,θ はガボール関数の模様の方向,ψ は位相オフセット,γ は空間ア
スペクト比を表す.ガボールフィルタの一例を Fig. 3.3.5 に示す.ガボールフィルタは画像処理でよ
く使用されており,指紋認証や虹彩認識などの生体認証,テクスチャのパターン認識などに用いられ
ている.また,2 次元ガボールフィルタにより一次視覚野の細胞の活動がモデル化できることが示さ
れている.
Fig. 3.3.4 Gabor 関数(自作)
Fig. 3.3.5 ガボールフィルタの一例(参考文献 4) より引用)
9
参考文献
1) ディジタル画像処理編集委員会. ディジタル画像処理. CG-ARTS 協会, 第 2 版, 2006.
2) 石田隆行, 大倉保彦, 青山正人, 川下郁生. C 言語で学ぶ医用画像処理. オーム社, 第 1 版, 2006.
3) 三谷政昭. やり直しのための通信数学:フーリエ変換からウェーブレット変換へ. CQ 出版社, 初
版, 2008.
4) ガボールフィルタ. http://ja.wikipedia.org/wiki/%E3%82%AC%E3%83%9C%E3%83%BC%E3%83%
AB%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF.
10