フーリエ変換¶
目的¶
- このチュートリアルでは
- OpenCVを使って画像のフーリエ変換を計算する方法を学びます．
- NumpyのFFTを使う方法を学びます．
- フーリエ変換を使ったアプリケーションについて学びます．
- 以下の関数の使い方を学びます : cv2.dft(), cv2.idft() etc
理論¶
フーリエ変換は種々のフィルタの周波数特性を解析するために使われます．画像に対しては 2次元離散フーリエ変換 (DFT) を使って周波数領域に変換します．高速化されたアルゴリズムである 高速フーリエ変換 (FFT) はDFTの計算に使います．これらのアルゴリズムの詳細については信号処理や画像処理の教科書を参照してください． 補足資料 の章に幾つか参考文献を挙げています．
正弦波を 
と書きます．ここで 
は信号の周波数を表します．もしこの信号を周波数領域で観測すると，周波数 の点にspikeが見られます．離散信号を形成するために信号を標本化すると，同じ周波数領域での信号を得られますが，周波数領域での信号は [- pi, pi] の範囲もしくは ![[0,2\pi]](../../../../_images/math/afd1cc9f6085210dc7aea0e591b90df796cff56c.png)
の範囲での周期性を持つ信号とみなされます(N点DFTであれば ![[0,N]](../../../../_images/math/2efe13038e0b745725e05c174d88172d3b57de17.png)
の範囲)．画像を2方向に標本化された信号とみなすことができます．横方向と縦方向にフーリエ変換を計算すれば，画像を周波数領域で表現できます．
より直観的に言うと，ある正弦波信号の振幅変化が短時間に早く起こると高周波な信号と言います．振幅変化が遅ければ低周波信号と呼びます．全く同じ考えを画像に対して拡張しましょう．画像中で振幅変化が急激に生じる場所はどこでしょうか?それはエッジやノイズです．つまり，エッジやノイズは画像の高周波成分に対応しています．振幅の変化が大きくなければ低周波成分になります( 周波数変換の直観的な説明をしている資料を 補足資料 の欄に幾つか挙げておきます)．
それではフーリエ変換の計算方法を見ていきましょう．
Numpyを使ったフーリエ変換¶
まず初めにNumpyを使ったフーリエ変換の計算方法を見ていきましょう．NumpyはFFTを計算するための関数 np.fft.fft2() を用意しています．この関数は複素数型の配列を出力します．第1引数は入力画像をグレースケール画像として与えます．第2引数は出力配列のサイズを指定しますが，オプションです．指定するサイズが入力画像のサイズより大きければ入力画像はFFTの計算をする前にゼロパディングされます．入力画像のサイズより小さければ入力画像を切り取ります．何も指定されなければ出力配列のサイズは入力画像のサイズと同じになります．
実際に計算をしてみると，周波数領域の原点(直流成分)が画像の左上の角に位置するようになります．直流成分を画像中心に移動させたければ，スペクトル全体を 
両方向のずらす必要があります．この移動には np.fft.fftshift() 関数を使います(これでより解析がしやすくなります)．一度フーリエ変換を計算すれば，スペクトルの大きさが分かります．
import cv2
import numpy as np
from matplotlib import pyplot as plt
img = cv2.imread('messi5.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()
結果は以下のようになります:
中心に白い領域が集中している事が分かります．これは画像が低周波成分をより多く含んでいることを意味します．
これでフーリエ変換を見れました．これでハイパスフィルタといった周波数領域での処理ができるようになります．低周波成分に対して矩形windowを使ったマスク処理をすることによって低周波成分を取り除く事が出来ます．それから np.fft.ifftshift() 関数を使って直流成分の位置を画像の左上に戻してから np.ifft2() 関数を使って逆フーリエ変換を適用します．最終的な結果は複素数型の配列になるので，その絶対値をとります．
rows, cols = img.shape
crow,ccol = rows/2 , cols/2
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)
plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])
plt.show()
結果は以下のようになります:
結果を見ると，ハイパスフィルタによって画像中のエッジが検出できていることが分かります．これは画像の勾配のチュートリアルで見た結果そのものですね．画像の大半の成分が低周波領域に存在する事も分かります．とにかく，これでNumpyを使って離散フーリエ変換（DFT)，逆離散フーリエ変換（IDFT)をする方法を学びました．次はOpenCVを使った方法を学びます．
結果を注意深く観察すると，特にJET色で描いた最後の画像を見ると，疑似輪郭(赤い矢印で示した部分です)が生じているのが分かります．波のような構造が見えますが，これを リンギング効果 と呼びます．マスク処理の際に矩形ウィンドウを使ったことが原因です．このマスクはフーリエ変換によってsinc関数になるためこのような結果になってしまいます．よって，矩形ウィンドウはフィルタリングには使いません．Gaussianウィンドウなどがより適しているでしょう．
OpenCVを使ったフーリエ変換¶
OpenCVはDFTを行う cv2.dft() とIDFTを行う cv2.idft() という関数を用意しています．Numpyと同様複素数型の配列を返しますが，2チャンネルの配列としてです．最初のチャンネルは結果の実部，二つ目のチャンネルが虚部に対応しています．入力画像は np.float32 型に変換される必要があります．それではどのように計算するのか見てみましょう．
import numpy as np
import cv2
from matplotlib import pyplot as plt
img = cv2.imread('messi5.jpg',0)
dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()
Note
cv2.cartToPolar() 関数を使い振幅と位相の両方を取得でいます．
次にするべきはIDFTです．前セッションでハイパスフィルタを試したので，ここではローパスフィルタ(高周波成分の除去)を試しましょう．ローパスフィルタは画像にボケを加えます．まず初めに低周波領域に高い値を持ち，高周波領域が0となるマスクを作成します．
rows, cols = img.shape
crow,ccol = rows/2 , cols/2
# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1
# apply mask and inverse DFT
fshift = dft_shift*mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()
結果を見てみましょう:
Note
いつも通り，OpenCVの関数である cv2.dft() と cv2.idft() はNumpyの処理より高速です．しかし，Numpyの関数の方が扱いやすいインターフェースを提供しています．パフォーマンスについては次の章を読んでください．
DFTのパフォーマンス最適化¶
DFTの計算に向いている配列サイズがあります．配列のサイズが2の累乗の時に最も高速に動きます．配列のサイズが2,3,5の積で表される時も効率的に計算されます．自分の実装に不安を感じるのであれば，DFTを適用する前に配列のサイズをゼロパディング等で最適なサイズに変更すると良いでしょう．OpenCVを使うのであればゼロパディングは手動で行う必要があります．一方で，NumpyであればFFTを計算する時の配列のサイズを指定すれば自動的にゼロパディングが行われます．
この最適な値をどのように計算すれば良いでしょうか?OpenCVが提供している cv2.getOptimalDFTSize() 関数を使います． cv2.dft() と np.fft.fft2() 関数のどちらにでも使えます．IPｙてょｎのmagic commandである %timeit を使ってパフォーマンスを調べてみましょう．
In [16]: img = cv2.imread('messi5.jpg',0)
In [17]: rows,cols = img.shape
In [18]: print rows,cols
342 548
In [19]: nrows = cv2.getOptimalDFTSize(rows)
In [20]: ncols = cv2.getOptimalDFTSize(cols)
In [21]: print nrows, ncols
360 576
配列のサイズが(342,548)から(360, 576)へ変更されたのが分かるでしょうか．ゼロパディングをしてからDFTをしてみましょう．ゼロパディングをするには，全要素の値を0にした新しい配列を用意し，元の画像データをコピーします．
nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img
もしくは cv2.copyMakeBorder() 関数を使ってもゼロパディングをできます:
right = ncols - cols
bottom = nrows - rows
bordertype = cv2.BORDER_CONSTANT #just to avoid line breakup in PDF file
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)
ようやくDFTを計算できます．Numpyの関数のパフォーマンスを比較してみましょう:
In [22]: %timeit fft1 = np.fft.fft2(img)
10 loops, best of 3: 40.9 ms per loop
In [23]: %timeit fft2 = np.fft.fft2(img,[nrows,ncols])
100 loops, best of 3: 10.4 ms per loop
4倍高速になりました．同様のパフォーマンスの計測をOpenCVでも行ってみましょう．
In [24]: %timeit dft1= cv2.dft(np.float32(img),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 13.5 ms per loop
In [27]: %timeit dft2= cv2.dft(np.float32(nimg),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 3.11 ms per loop
こちらも4倍高速化されました．OpenCVの関数を使うとNumpyの関数を使うより3倍ほど速く計算が出来ているのが分かります．逆フーリエ変換でも同じテストをしてみましょう．これはあなたの課題として残しておきます．
何故Laplacianがハイパスフィルタなのか?¶
フォーラムでも似たような質問が聞かれています．質問とは，何故Laplacianフィルタがハイパスフィルタなのか?何故Sobelフィルタがハイパスフィルタなのか?です．最初の答えはフーリエ変換の観点から説明します．Laplacianフィルタのフーリエ変換を計算し，解析してみましょう:
import cv2
import numpy as np
from matplotlib import pyplot as plt
# simple averaging filter without scaling parameter
mean_filter = np.ones((3,3))
# creating a guassian filter
x = cv2.getGaussianKernel(5,10)
gaussian = x*x.T
# different edge detecting filters
# scharr in x-direction
scharr = np.array([[-3, 0, 3],
[-10,0,10],
[-3, 0, 3]])
# sobel in x direction
sobel_x= np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]])
# sobel in y direction
sobel_y= np.array([[-1,-2,-1],
[0, 0, 0],
[1, 2, 1]])
# laplacian
laplacian=np.array([[0, 1, 0],
[1,-4, 1],
[0, 1, 0]])
filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
filter_name = ['mean_filter', 'gaussian','laplacian', 'sobel_x', \
'sobel_y', 'scharr_x']
fft_filters = [np.fft.fft2(x) for x in filters]
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
mag_spectrum = [np.log(np.abs(z)+1) for z in fft_shift]
for i in xrange(6):
plt.subplot(2,3,i+1),plt.imshow(mag_spectrum[i],cmap = 'gray')
plt.title(filter_name[i]), plt.xticks([]), plt.yticks([])
plt.show()
結果を見てください:
結果画像を見ると，各カーネルがどの周波数帯域をブロックするのか，どの周波数帯域を通過させるのか分かります．ここから上記のデモで試したカーネルはハイパスフィルタもしくはローパスフィルタのどちらかになることが分かります．
補足資料¶
- フーリエ変換の理論の直観的な説明(英語) by Steven Lehar
- フーリエ変換(英語) at HIPR
- 画像における周波数領域が何を意味するのか(英語)
