Pythonで画像処理（2次元FFTとパワースペクトル）

画像のパワースペクトル（2次元FFTの絶対値の2乗）を画像で出力するプログラムをPythonで書いた。

とにかく、コードを載せる。

spectrum.py

import sys
from PIL import Image
import numpy as np

if len(sys.argv) != 3:
    print('Usage: python spectrum.py input output')
    sys.exit(1)

img_in = Image.open(sys.argv[1])
x = np.array(img_in.convert('L'))

X = np.fft.fft2(x)
X = np.fft.fftshift(X)

X_abs = np.absolute(X)
X_abs[X_abs < 1] = 1  # Assign 1 to elements (<1)
P = np.log10(X_abs)  # Power spectrum

# Normalize to [0, 255]
P_norm = P/np.amax(P)
y = np.uint8(np.around(P_norm*255))

img_out = Image.fromarray(y)
img_out.save(sys.argv[2])

まず、convert()でグレースケールに変換。特定の色成分だけ取り出して見ることもあると思うが、ここでは（カラー画像でも）輝度値で見ている。np.fft.fft2()が2次元FFTで、np.fft.fftshift()は何かというと、すっかり忘れていたのだが以前C言語で実装したときの説明があったのでこれでいいや。

yaritakunai.hatenablog.com

2次元FFTした結果は、直流成分が配列の左上にあります。
パワースペクトルを見やすくするために、まずnp.fft.fftshift()を使って直流成分を配列の中心に移動させます。
np.fft.fftshift()は、配列の第1象限と第4象限、第2象限と第3象限をそれぞれ入れ替えています。

（一部改変）

その後の処理はちと違う。輝度は非負値のため、そのFFTは直流成分に大きな値を持つ。パワースペクトルなら2乗しているのでなおさらである。そのため、線形スケールだと高々256段階ではほとんどの成分が黒で潰れてしまう。そこで、ログスケールにしている。

パワーの定義式では係数に20が掛かるが、どうせ正規化するので計算しない。また、対数は1以下の引数でマイナスになるので、1（つまり、ログとって0）で潰している。こういうとき、NumPyのインデックス記法は便利。

多少やっつけ記事だが、出力結果を貼って終わりとする。（左は入力画像）

実践 コンピュータビジョン