Pythonで画像をフーリエ変換してフィルタをかけてみた | 診療放射線技師がPythonをはじめました。

画像をフーリエ変換して何をするかといえば、やっぱりフィルタですよね。

ということで今回は画像をフーリエ変換してローパスフィルタ、ハイパスフィルタを試してみたのでそれを記事にします。

今回使用する画像は以下の画像です。

Contents

1 ローパスフィルタ
2 ハイパスフィルタ
3 最後に

ローパスフィルタ

まず初めに単純に高周波成分をカットして低周波成分だけを通過させるローパスフィルタを試してみました。

やり方は簡単で、フーリエ変換して位相シフトした画像と同じサイズの配列を作成し低周波成分にあたる部分を１、高周波成分にあたる部分を０とします。

その配列と、フーリエ変換後の配列を掛け算することでできます。

import numpy as np
import matplotlib.pyplot as plt

img = pydicom.dcmread(‘CT000034′)
pix_arr = img.pixel_array           #　変数pix_arrにピクセルデータを格納

fft_arr = np.fft.fft2(pix_arr)      # 変数fft_arrにpix_arrをフーリエ変換したデータを格納
fft_arr_shift = np.fft.fftshift(fft_arr)

filter = np.zeros(fft_arr.shape)
center_x = fft_arr.shape&#91;1]//2
center_y = fft_arr.shape&#91;0]//2

filter&#91;center_y-50: center_y+50 ,center_x-50: center_x+50 ] = 1

fft_arr_shift_filter = fft_arr_shift * filter

fft_filter = np.fft.fftshift(fft_arr_shift_filter)   # 位相シフトを戻す
fft_filter = np.fft.ifft2(fft_filter)                 #逆フーリエ変換

fig = plt.figure()
ax1 = fig.add_subplot(2,3,1)
ax2 = fig.add_subplot(2,3,2)
ax3 = fig.add_subplot(2,3,3)
ax4 = fig.add_subplot(2,3,4)
ax5 = fig.add_subplot(2,3,4)
ax6 = fig.add_subplot(2,3,6)

ax1.imshow(np.abs(fft_arr_shift.real), cmap=’gray’, vmax =500000)
ax2.imshow(np.abs(filter), cmap=’gray’)
ax3.imshow(np.abs(fft_arr_shift_filter), cmap=’gray’, vmax =500000)
ax4.imshow(pix_arr, cmap=’gray’)
ax6.imshow(np.abs(fft_filter.real), cmap=’gray’)

plt.show()

import numpy as np

import matplotlib.pyplot as plt

img = pydicom.dcmread(‘CT000034’)

pix_arr = img.pixel_array #　変数pix_arrにピクセルデータを格納

fft_arr = np.fft.fft2(pix_arr) # 変数fft_arrにpix_arrをフーリエ変換したデータを格納

fft_arr_shift = np.fft.fftshift(fft_arr)

filter = np.zeros(fft_arr.shape)

center_x = fft_arr.shape[1]//2

center_y = fft_arr.shape[0]//2

filter[center_y-50: center_y+50 ,center_x-50: center_x+50 ] = 1

fft_arr_shift_filter = fft_arr_shift * filter

fft_filter = np.fft.fftshift(fft_arr_shift_filter) # 位相シフトを戻す

fft_filter = np.fft.ifft2(fft_filter) #逆フーリエ変換

fig = plt.figure()

ax1 = fig.add_subplot(2,3,1)

ax2 = fig.add_subplot(2,3,2)

ax3 = fig.add_subplot(2,3,3)

ax4 = fig.add_subplot(2,3,4)

ax5 = fig.add_subplot(2,3,4)

ax6 = fig.add_subplot(2,3,6)

ax1.imshow(np.abs(fft_arr_shift.real), cmap=‘gray’, vmax =500000)

ax2.imshow(np.abs(filter), cmap=‘gray’)

ax3.imshow(np.abs(fft_arr_shift_filter), cmap=‘gray’, vmax =500000)

ax4.imshow(pix_arr, cmap=‘gray’)

ax6.imshow(np.abs(fft_filter.real), cmap=‘gray’)

plt.show()

画像が小さくてわかりずらいかもしれませんが、画像のシャープさがなくなっているのが分かります。（クリックすると拡大します。）

フィルタの１の領域をもう少し小さくしてみます。上記コードの14行目を書き換えます。

filter[center_y-5: center_y+5 ,center_x-5: center_x+5 ] = 1

1	filter[center_y–5: center_y+5 ,center_x–5: center_x+5 ] = 1

前のフィルタの1/10にしたらアイソトープのような分解能の画像になってしまいました。。。。。

ハイパスフィルタ

今度はハイパスフィルタをやってみます。

先ほどは、0で初期値設定した配列を作成しましたが、今度は１で初期化した配列を作成し、中心部を0で埋めた配列を作成してハイパスフィルタを作成します。

import numpy as np
import matplotlib.pyplot as plt

img = pydicom.dcmread(‘CT000033′)
pix_arr = img.pixel_array           #　変数pix_arrにピクセルデータを格納

fft_arr = np.fft.fft2(pix_arr)      # 変数fft_arrにpix_arrをフーリエ変換したデータを格納
fft_arr_shift = np.fft.fftshift(fft_arr)

filter = np.ones(fft_arr.shape)
center_x = fft_arr.shape&#91;1]//2
center_y = fft_arr.shape&#91;0]//2

filter&#91;center_y-100: center_y+100 ,center_x-100: center_x+100 ] = 0

fft_arr_shift_filter = fft_arr_shift * filter

fft_filter = np.fft.fftshift(fft_arr_shift_filter)   # 位相シフトを戻す
fft_filter = np.fft.ifft2(fft_filter)                 #逆フーリエ変換

fig = plt.figure()
ax1 = fig.add_subplot(2,3,1)
ax2 = fig.add_subplot(2,3,2)
ax3 = fig.add_subplot(2,3,3)
ax4 = fig.add_subplot(2,3,4)
ax5 = fig.add_subplot(2,3,4)
ax6 = fig.add_subplot(2,3,6)

ax1.imshow(np.abs(fft_arr_shift.real), cmap=’gray’, vmax =500000)
ax2.imshow(np.abs(filter), cmap=’gray’)
ax3.imshow(np.abs(fft_arr_shift_filter), cmap=’gray’, vmax =500000)
ax4.imshow(pix_arr, cmap=’gray’)
ax6.imshow(np.abs(fft_filter.real), cmap=’gray’)

plt.show()