フーリエ変換は音声データに対して用いられることが多い手法だけど、画像データにも応用が効く。 音声データの場合、フーリエ変換を使うことで時間領域の情報を周波数領域の情報に直せる。 それに対し、画像データでは空間領域の情報を周波数領域の情報に直すことになる。 つまり、画像データの濃淡が複数の波形の合成によって作られていると見なす。 今回は、画像データをフーリエ変換して、周波数領域の情報にフィルタをかけたり元に戻したりして遊んでみる。
使った環境は次のとおり。
$ sw_vers ProductName: Mac OS X ProductVersion: 10.14.6 BuildVersion: 18G5033 $ python -V Python 3.7.7
もくじ
- もくじ
- 下準備
- 周波数領域の情報を可視化してみる
- ローパスフィルタ (Low Pass Filter) をかけて復元してみる
- ハイパスフィルタ (High Pass Filter) をかけて復元してみる
- カラー画像を扱う場合
下準備
あらかじめ、必要なパッケージをインストールしておく。
$ pip install pillow matplotlib
周波数領域の情報を可視化してみる
最初に、画像データをフーリエ変換して周波数領域で可視化してみる。
画像データは、二次元の NumPy 配列として読み込んだ上で np.fft.fft2() に渡すことでフーリエ変換できる。
また、周波数領域のデータは np.fft.ifft2() を使うことで空間領域のデータに戻せる。
以下のサンプルコードでは、元の画像と、周波数領域でのパワースペクトル、そして逆変換することで元に戻した画像を可視化している。 なお、読み込む画像は適当に用意しよう。
#!/usr/bin/env python3 # -*- coding: utf-8 -*- from PIL import Image import numpy as np from matplotlib import pyplot as plt def main(): # 画像を読み込む img = Image.open('lena.png') # グレイスケールに変換する gray_img = img.convert('L') # NumPy 配列にする f_xy = np.asarray(gray_img) # 2 次元高速フーリエ変換で周波数領域の情報を取り出す f_uv = np.fft.fft2(f_xy) # 画像の中心に低周波数の成分がくるように並べかえる shifted_f_uv = np.fft.fftshift(f_uv) # パワースペクトルに変換する magnitude_spectrum2d = 20 * np.log(np.absolute(shifted_f_uv)) # 元の並びに直す unshifted_f_uv = np.fft.fftshift(shifted_f_uv) # 2 次元逆高速フーリエ変換で空間領域の情報に戻す i_f_xy = np.fft.ifft2(unshifted_f_uv).real # 実数部だけ使う # 上記を画像として可視化する fig, axes = plt.subplots(1, 3, figsize=(8, 4)) # 枠線と目盛りを消す for ax in axes: for spine in ax.spines.values(): spine.set_visible(False) ax.set_xticks([]) ax.set_yticks([]) # 元画像 axes[0].imshow(f_xy, cmap='gray') axes[0].set_title('Input Image') # 周波数領域のパワースペクトル axes[1].imshow(magnitude_spectrum2d, cmap='gray') axes[1].set_title('Magnitude Spectrum') # FFT -> IFFT した画像 axes[2].imshow(i_f_xy, cmap='gray') axes[2].set_title('Reversed Image') # グラフを表示する plt.show() if __name__ == '__main__': main()
上記をファイルに保存して実行する。
$ python imgfft2d.py
すると、次のようなグラフが得られる。
真ん中のプロットが、元となった画像データの周波数領域での表現になる。 この表現では、中心に近いほど低い周波数・遠いほど高い周波数の成分を含む。 白い部分ほど成分が多いことを示しているため、この画像は低周波の成分が比較的多いように見える。 また、右の画像を見ると、ちゃんと周波数領域のデータから空間領域のデータに復元できていることがわかる。
ローパスフィルタ (Low Pass Filter) をかけて復元してみる
先ほど述べた通り、周波数領域の表現では中心に近いほど低い周波数・遠いほど高い周波数の成分を含んでいる。 つまり、中心部分だけを取り出すような演算をすると、画像を構成する波形にローパスフィルタをかけることができる。
以下のサンプルコードでは、周波数領域のデータに対して中心部分を取り出すフィルタをかけることで低周波成分だけを抽出した。 そして、フィルタしたデータを逆フーリエ変換することで画像に戻している。
#!/usr/bin/env python3 # -*- coding: utf-8 -*- from PIL import Image from PIL import ImageDraw import numpy as np from matplotlib import pyplot as plt def main(): # 画像を読み込む img = Image.open('lena.png') # グレイスケールに変換する gray_img = img.convert('L') # NumPy 配列にする f_xy = np.asarray(gray_img) # 2 次元高速フーリエ変換で周波数領域の情報を取り出す f_uv = np.fft.fft2(f_xy) # 画像の中心に低周波数の成分がくるように並べかえる shifted_f_uv = np.fft.fftshift(f_uv) # フィルタ (ローパス) を用意する x_pass_filter = Image.new(mode='L', # 8-bit pixels, black and white size=(shifted_f_uv.shape[0], shifted_f_uv.shape[1]), color=0, # default black ) # 中心に円を描く draw = ImageDraw.Draw(x_pass_filter) # 円の半径 ellipse_r = 50 # 画像の中心 center = (shifted_f_uv.shape[0] // 2, shifted_f_uv.shape[1] // 2) # 円の座標 ellipse_pos = (center[0] - ellipse_r, center[1] - ellipse_r, center[0] + ellipse_r, center[1] + ellipse_r) draw.ellipse(ellipse_pos, fill=255) # フィルタ filter_array = np.asarray(x_pass_filter) # フィルタを適用する filtered_f_uv = np.multiply(shifted_f_uv, filter_array) # パワースペクトルに変換する magnitude_spectrum2d = 20 * np.log(np.absolute(filtered_f_uv)) # 元の並びに直す unshifted_f_uv = np.fft.fftshift(filtered_f_uv) # 2 次元逆高速フーリエ変換で空間領域の情報に戻す i_f_xy = np.fft.ifft2(unshifted_f_uv).real # 実数部だけ使う # 上記を画像として可視化する fig, axes = plt.subplots(1, 4, figsize=(12, 4)) # 枠線と目盛りを消す for ax in axes: for spine in ax.spines.values(): spine.set_visible(False) ax.set_xticks([]) ax.set_yticks([]) # 元画像 axes[0].imshow(f_xy, cmap='gray') axes[0].set_title('Input Image') # フィルタ画像 axes[1].imshow(filter_array, cmap='gray') axes[1].set_title('Filter Image') # フィルタされた周波数領域のパワースペクトル axes[2].imshow(magnitude_spectrum2d, cmap='gray') axes[2].set_title('Filtered Magnitude Spectrum') # FFT -> Band-pass Filter -> IFFT した画像 axes[3].imshow(i_f_xy, cmap='gray') axes[3].set_title('Reversed Image') # グラフを表示する plt.show() if __name__ == '__main__': main()
上記を実行してみよう。
$ python imglpf.py
すると、次のようなグラフが得られる。 分かりにくいかもしれないけど、復元した画像はちょっとぼやけている。 また、ぼやけ方も縞模様になっていて波を感じるものとなった。 これは、画像を構成するすべての周波数の中から、高周波の成分がフィルタによって取り除かれたことで生じている。
中心のごく一部だけを取り出している (つまり、多くの情報が失われている) のに、それっぽい画像になるのはなんとも面白い。 なお、フィルタに使う円の半径を小さくすれば、それだけ高周波の成分が少なくなって、より分かりやすくぼやける。
ハイパスフィルタ (High Pass Filter) をかけて復元してみる
次は、同様に高周波成分だけを取り出すハイパスフィルタをかけてみよう。 つまり、中心だけを取り除くことになる。 先ほどとのコードの違いは、適用するフィルタの周波数特性が違うだけ。
#!/usr/bin/env python3 # -*- coding: utf-8 -*- from PIL import Image from PIL import ImageDraw import numpy as np from matplotlib import pyplot as plt def main(): # 画像を読み込む img = Image.open('lena.png') # グレイスケールに変換する gray_img = img.convert('L') # NumPy 配列にする f_xy = np.asarray(gray_img) # 2 次元高速フーリエ変換で周波数領域の情報を取り出す f_uv = np.fft.fft2(f_xy) # 画像の中心に低周波数の成分がくるように並べかえる shifted_f_uv = np.fft.fftshift(f_uv) # フィルタ (ハイパス) を用意する x_pass_filter = Image.new(mode='L', # 8-bit pixels, black and white size=(shifted_f_uv.shape[0], shifted_f_uv.shape[1]), color=255, # default white ) # 中心に円を描く draw = ImageDraw.Draw(x_pass_filter) # 円の半径 ellipse_r = 50 # 画像の中心 center = (shifted_f_uv.shape[0] // 2, shifted_f_uv.shape[1] // 2) # 円の座標 ellipse_pos = (center[0] - ellipse_r, center[1] - ellipse_r, center[0] + ellipse_r, center[1] + ellipse_r) draw.ellipse(ellipse_pos, fill=0) # フィルタ filter_array = np.asarray(x_pass_filter) # フィルタを適用する filtered_f_uv = np.multiply(shifted_f_uv, filter_array) # パワースペクトルに変換する magnitude_spectrum2d = 20 * np.log(np.absolute(filtered_f_uv)) # 元の並びに直す unshifted_f_uv = np.fft.fftshift(filtered_f_uv) # 2 次元逆高速フーリエ変換で空間領域の情報に戻す i_f_xy = np.fft.ifft2(unshifted_f_uv).real # 実数部だけ使う # 上記を画像として可視化する fig, axes = plt.subplots(1, 4, figsize=(12, 4)) # 枠線と目盛りを消す for ax in axes: for spine in ax.spines.values(): spine.set_visible(False) ax.set_xticks([]) ax.set_yticks([]) # 元画像 axes[0].imshow(f_xy, cmap='gray') axes[0].set_title('Input Image') # フィルタ画像 axes[1].imshow(filter_array, cmap='gray') axes[1].set_title('Filter Image') # フィルタされた周波数領域のパワースペクトル axes[2].imshow(magnitude_spectrum2d, cmap='gray') axes[2].set_title('Filtered Magnitude Spectrum') # FFT -> Band-pass Filter -> IFFT した画像 axes[3].imshow(i_f_xy, cmap='gray') axes[3].set_title('Reversed Image') # グラフを表示する plt.show() if __name__ == '__main__': main()
上記を保存して実行する。
$ python imghpf.py
すると、今度は次のようなグラフが得られる。 またもや分かりにくいけど、復元した画像は元の画像の輪郭だけがぼんやりと浮かび上がったものとなっている。 これが画像を構成する高周波の成分で、先ほどの復元画像から取り除かれたものと解釈できる。
カラー画像を扱う場合
なお、カラー画像を扱う場合、それぞれのチャネルごとに処理する必要がある。 以下のサンプルコードでは、最初のサンプルコードを RGB のチャネルでそれぞれ適用している。
#!/usr/bin/env python3 # -*- coding: utf-8 -*- from itertools import chain from PIL import Image import numpy as np from matplotlib import pyplot as plt def main(): # 画像を読み込む img = Image.open('lena.png') # NumPy 配列にする f_xy_rgb = np.asarray(img) # 画像として可視化する fig, axes = plt.subplots(3, 5, figsize=(12, 6)) # 枠線と目盛りを消す for ax in chain.from_iterable(axes): for spine in ax.spines.values(): spine.set_visible(False) ax.set_xticks([]) ax.set_yticks([]) # 元画像 axes[1][0].imshow(f_xy_rgb) axes[1][0].set_title('Input Image (RGB)') # 逆変換した画像を入れる領域を用意しておく i_f_xy_rgb = np.empty_like(f_xy_rgb) # RGB ごとに処理する channel_names = ('R', 'G', 'B') channel_indices = range(f_xy_rgb.shape[2]) zipped_channel_info = zip(channel_names, channel_indices) for channel_name, channel_index in zipped_channel_info: # 各チャネルごとに配列として取り出す f_xy = f_xy_rgb[:, :, channel_index] # 各チャネルの元画像 axes[channel_index][1].imshow(f_xy, cmap='gray') axes[channel_index][1].set_title(f'Input Image ({channel_name})') # 2 次元高速フーリエ変換で周波数領域の情報を取り出す f_uv = np.fft.fft2(f_xy) # 画像の中心に低周波数の成分がくるように並べかえる shifted_f_uv = np.fft.fftshift(f_uv) # パワースペクトルに変換する magnitude_spectrum2d = 20 * np.log(np.absolute(shifted_f_uv)) # 元の並びに直す unshifted_f_uv = np.fft.fftshift(shifted_f_uv) # 2 次元逆高速フーリエ変換で空間領域の情報に戻す i_f_xy = np.fft.ifft2(unshifted_f_uv).real # 実数部だけ使う # 周波数領域のパワースペクトル axes[channel_index][2].imshow(magnitude_spectrum2d, cmap='gray') axes[channel_index][2].set_title(f'Magnitude Spectrum ({channel_name})') # FFT -> Band-pass Filter -> IFFT した画像 axes[channel_index][3].imshow(i_f_xy, cmap='gray') axes[channel_index][3].set_title(f'Reversed Image ({channel_name})') # 逆変換したチャネルを保存しておく i_f_xy_rgb[:, :, channel_index] = i_f_xy # 逆変換した RGB 画像 axes[1][4].imshow(i_f_xy_rgb) axes[1][4].set_title('Reversed Image (RGB)') # グラフを表示する plt.show() if __name__ == '__main__': main()
上記を保存して実行する。
$ python imgfftc.py
すると、次のようなグラフが得られる。
ちゃんと RGB それぞれのチャネルの情報から、カラー画像が復元できていることがわかる。
めでたしめでたし。
![ディジタル画像処理[改訂第二版]](https://m.media-amazon.com/images/I/41FCF736uWL._SL500_.jpg "ディジタル画像処理[改訂第二版]")
