【やってみた】Pythonでフーリエ変換を使って画像処理してみた

株式会社神戸デジタル・ラボ DataIntelligenceチームの高木です。

今回は、画像の空間と異なる空間上で処理を行うフーリエ変換を利用した画像処理について解説していきたいと思います。

具体例を通して理解するため、本記事ではこんなテーマを用意しました。

「檻に閉じ込められたデジごんを救出！」

早速、内容に入っていきましょう！

デジごんとは？
デジごんが檻に閉じ込められた？
画像におけるフーリエ変換
檻から救出大作戦
まとめ

デジごんとは？

おそらく皆さんがテーマを見て最初に思われたのは「デジごんって何だ？」ですよね！

こちらは弊社の公式キャラクターになります。

のほほんとした、顔がなんとも言えません。

『デジごんのLINEスタンプ販売中！』

store.line.me

カラーのデジごんも捨てがたいところですが、今回は実装コードの簡略化の都合上デジごんには白黒になってもらいまいしょう！！ついでに、ポージングもいけてる感じのものを採用します！

デジごんが檻に閉じ込められた？

下の画像をご覧ください。どうやら、デジごんが檻に閉じ込められたようです。。。

画像から檻を取り除いて、閉じ込められたデジごんを解放してあげないと！！！

そんな時に有効なのが画像におけるフーリエ変換です。

画像におけるフーリエ変換

この章では、フーリエ変換における数学的な部分の説明を少ししていきます。数式の理解ができなくても、以降の章には影響しませんので、数式が苦手な方は流し読みで大丈夫です！

画像をある関数 $f(x, y)$ とするとき、フーリエ変換 $F(u, v)$ は以下のように定義されます。

$F(u, v) = \iint_{-\infty}^{ \infty} f(x, y) \hspace{2mm} \text{exp} \lbrace-i2\pi(ux + vy)\rbrace dxdy$

（ $i = \sqrt{-1}$ , 　 $u$ ： $x$ 方向の空間周波数,　 $v$ ： $y$ 方向の空間周波数）

数式で見るとなんだか難しそうなんですが、フーリエ変換というのは $x, y$ で表される画像を表す関数 $f(x, y)$ が、 $u, v$ で表される別の空間上で関数 $F(u, v)$ といった別の形で表現できるということを表します。（ $x, y$ で表される空間は空間領域、 $u, v$ で表される空間は周波数領域と呼ばれます）

ちなみに、空間周波数とは空間的な構造における繰り返しの数を指します。感覚的にとらえることはなかなか難しいところでもあるので、そういう値なんだなあという程度の認識で大丈夫です！

参考）周波数は、下図のような周期の繰り返しにおいて1秒間における周期の数のことを示します。

また、関数 $F(u, v)$ は以下のようなフーリエ逆変換によって、元の画像を表す関数 $f(x, y)$ に戻すことができます。

$f(x, y) = \iint_{-\infty}^{ \infty} F(u, v) \hspace{2mm} \text{exp} \lbrace i2\pi(ux + vy)\rbrace dudv$

これらの関係をまとめるとこんな感じですかね！

周波数領域に関連する値にはそれぞれ名前がついています。

振幅スペクトル　 $|F(u, v)| = \sqrt{Re\lbrace F(u, v) \rbrace^{2} + Im\lbrace F(u, v) \rbrace^{2}}$

位相スペクトル　 $arg \lbrace F(u, v) \rbrace = tan^{-1} \dfrac{Im\lbrace F(u, v) \rbrace}{Re\lbrace F(u, v) \rbrace}$

パワースペクトル　 $|F(u, v)|^{2}$

檻から救出大作戦

まずは、今回のプログラムを実行するために必要なライブラリをインポートしておきましょう！

import numpy as np
from matplotlib import pyplot as plt
import cv2

檻にいるデジごんをフーリエ変換してみよう

檻に閉じ込められたデジごんの画像に対して、先ほど紹介したフーリエ変換をしてみます。

数式で見ると難しそうだったフーリエ変換もプログラムではとても簡単です。np.fft.fft2の部分がフーリエ変換に当たり、それ以外の部分は画像の読み込みとフーリエ変換ができるような形式（numpy形式）に変更しているだけになります！

main_filename = "degigon.png"
main_image = cv2.imread(main_filename, cv2.IMREAD_GRAYSCALE)
main_image = np.asarray(main_image)
main_fq = np.fft.fft2(main_image) # フーリエ変換

上記のコードでフーリエ変換は完了したのですが、変換後がどのようになっているのかは気になる部分ですよね？そのためにパワースペクトル画像というものを作成します！前の章の最後で登場したパワースペクトルの値を可視化した画像ですね。画像としてみてあげることで、値の大小を色の濃淡として捉えることができます。

main_fq_shifted = np.fft.fftshift(main_fq) # 位置調整
main_spectrum = 20 * np.log(np.abs(main_fq_shifted)) # 数値調整
main_spectrum[main_spectrum == -np.inf] = np.nanmin(main_spectrum[main_spectrum != -np.inf]) # 無限大に対する対応

上記のコードで大切な部分は、1行目です。np.fft.fftshift では低周波成分を中央に・高周波成分を四隅に移動するといったことをおこなっています。これを行う理由としては、視覚的にわかりやすい画像になるからですね。

2行目では、スペクトルの特徴を掴みやすいためにスケールの変換を行っています。また、3行目では2行目の計算の際に、マイナス無限大になってしまう値を最小値で代替する処理を行なっています。

では、実際に完成したスペクトル画像を見てみましょう！

plt.imshow(main_spectrum, cmap='gray')
plt.axis('off')
plt.show()

なんだこれはという感じですが、周波数領域を表すスペクトル画像ではデジごんがこんな風に見えてるという捉え方でOKです！

以上のコードをまとめると、以下のようになります。

main_image = cv2.imread(main_filename, cv2.IMREAD_GRAYSCALE)
main_image = np.asarray(main_image)
main_fq = np.fft.fft2(main_image) 

main_fq_shifted = np.fft.fftshift(main_fq)
main_spectrum = 20 * np.log(np.abs(main_fq_shifted))
main_spectrum[main_spectrum == -np.inf] = np.nanmin(main_spectrum[main_spectrum != -np.inf])

plt.imshow(main_spectrum, cmap='gray')
plt.axis('off')
plt.show()

檻だけの画像をフーリエ変換してみよう

デジごんが閉じ込められている檻だけの画像も先ほどと同様にスペクトル変換して、パワースペクトル画像を出力してみます！

sub_filename = "ori.png"
sub_image = cv2.imread(sub_filename, cv2.IMREAD_GRAYSCALE)
sub_image = np.asarray(sub_image)
sub_fq = np.fft.fft2(sub_image) 

sub_fq_shifted = np.fft.fftshift(sub_fq)
sub_spectrum = 20 * np.log(np.abs(sub_fq_shifted))
sub_spectrum[sub_spectrum == -np.inf] = np.nanmin(sub_spectrum[sub_spectrum != -np.inf])

plt.imshow(sub_spectrum, cmap='gray')
plt.axis('off')
plt.show()

檻の画像は、どうやらパワースペクトル画像においても檻っぽくなるようですね！

変換後の差分をとってフーリエ逆変換しよう

ここでは、デジごんが閉じ込められた画像におけるフーリエ変換後の値と、檻のみの画像におけるフーリエ変換後の値の差分を取ることによって、檻の情報を取り除くことを考えます。

プログラム的には、先ほど登場した変数のmain_fqとsub_fqの引き算を取るだけなのでとても簡単です！この計算結果のdiff_fqを、フーリエ逆変換します。そうすることで、周波数領域から画像を表す空間領域に戻すことができるのです！

diff_fq = main_fq - sub_fq
image_tm_inverse = np.fft.ifft2(diff_fq).real

あとは、これを可視化するだけ！ということで、早速結果をみてみましょう！

plt.imshow(image_tm_inverse, cmap='gray')
plt.axis('off')
plt.savefig('final.png')
plt.show()

どうでしょうか？デジごんと被っていた檻の箇所が一部消えてしまいましたが、無事に檻から解放されました！

まとめ

今回は、フーリエ変換を利用した画像処理を理解するために、檻に閉じ込められたデジごんの救出を例に取り上げました！

今回行った檻を取り除く処理は、フーリエ変換前の空間領域だけで行うことはとても難しいです。しかし、フーリエ変換後の周波数領域に空間を移動することで引き算というとても簡単な手法で解決することができました！

【参考】
元画像と檻の画像のピクセル同士の引き算で檻から出せることは出せるのですが、そのためには画像のどの位置に檻があるのかをあらかじめ知っておく必要があります。

このように、フーリエ変換をした空間では特定の周波数成分（檻の部分に値する部分）を取り除くことが比較的に容易にできます。

みなさんも色々な画像で、フーリエ変換後の空間での処理（周波数領域）を試してみてください！きっと、新たな発見があるはずです！

高木裕仁

データインテリジェンスチーム所属
データサイエンティスト。自然言語処理を中心としながら、その他の非構造化データや構造化データに関しても偏りなく扱います。こちらのブログでは、自然言語処理に関するトピックやAzureを中心としたクラウドを利用したデータ活用に関してのトピックを中心に様々な記事を発信していきます。