フーリエ変換をして、画像の特定の周波数成分を除去するプログラムの紹介（Pythonでの実装）

本日は画像のフーリエ変換を用いて画像を周波数成分に変換し、特定の周波数成分をカットしてから元に戻すことで、画像を変化させてみます。

画像のフーリエ変換

フーリエ変換は、ある時間領域の信号を周波数領域の信号に変換する変換手法です。

画像に対しては2次元離散フーリエ変換が使われますが、これを高速に計算する手段としてFFT（Fast Fourier Transform：高速フーリエ変換）が広く用いられています。

一般的な画像は、低周波成分が多く、高周波成分が少ない傾向にあります。

よって低周波成分だけを残しても、画像の多くの特徴は残ります。よって、低周波成分だけを残して保存することで、画像のデータサイズを圧縮することができます。

また、一般にホワイトノイズと呼ばれるようなノイズは、低周波・高周波によらず全周波数帯に一様にノイズ成分が載る特徴があります。

よって、高周波成分を積極的に削減することによって、ノイズの成分のみを多く削減することができます。

このように、信号の見方を変えて周波数成分で扱うことは、画像処理において非常に強力なツールとなり得ます。

画像の特定の周波数をカットするプログラム（Pythonでの実装）

以下に画像を入力にフーリエ変換を施し、特定の周波数を除去するプログラムを紹介します。今回はグレースケール画像に対してフーリエ変換を実施します。

• 動作環境：OpenCV 4.5.5

	import cv2
	import sys
	import math
	import numpy as np
	from matplotlib import pyplot as plt
	def main():
	print("OpenCV Version: " + str(cv2.__version__) + "\n")
	# Loading noisy image data
	filename = "noisy_image.png"
	image = cv2.imread(filename, cv2.IMREAD_GRAYSCALE)
	# Loading groundtruth image data
	filename_gt = "gt_image.png"
	image_gt = cv2.imread(filename_gt, cv2.IMREAD_GRAYSCALE)
	if image is None:
	print("Cannot find image : " + filename)
	sys.exit()
	if image_gt is None:
	print("Cannot find image : " + filename_gt)
	sys.exit()
	image_tm = np.asarray(image)
	image_tm_gt = np.asarray(image_gt)
	# FFT (time domain -> frequency domain)
	image_fq = np.fft.fft2(image_tm);
	image_fq_shifted = np.fft.fftshift(image_fq)
	# Frequency mask creation
	CUTOFF_FREQ = 50;
	size = image_tm.shape
	mask = np.zeros(size)
	length_x = size[1]
	length_y = size[0]
	center_x = size[1]/2
	center_y = size[0]/2
	for x in range(0,length_x):
	for y in range(0,length_y):
	if abs(x - center_x) < CUTOFF_FREQ and abs(y - center_y) < CUTOFF_FREQ:
	mask[x,y]=1
	else:
	mask[x,y]=1e-10
	# Frequency mask
	image_fq_shifted = image_fq_shifted * mask
	image_fq = np.fft.fftshift(image_fq_shifted)
	# IFFT (frequency domain -> time domain)
	image_tm_inverse = np.fft.ifft2(image_fq).real;
	# Calculating power spectrum
	magnitude_spectrum = 20 * np.log(np.abs(image_fq_shifted));
	# Evaluation with PSNR
	PSNREval(image_tm_gt, image_tm, 255)
	PSNREval(image_tm_gt, image_tm_inverse, 255)
	plt.subplot(221),plt.imshow(image_tm, cmap = 'gray')
	plt.title('Image (input)'), plt.xticks([]), plt.yticks([])
	plt.subplot(222),plt.imshow(mask, cmap = 'gray')
	plt.title('Frequency filter'), plt.xticks([]), plt.yticks([])
	plt.subplot(223),plt.imshow(image_tm_inverse, cmap = 'gray')
	plt.title('Image (inverse)'), plt.xticks([]), plt.yticks([])
	plt.subplot(224),plt.imshow(magnitude_spectrum, cmap = 'gray')
	plt.title('Magnitude spectrum'), plt.xticks([]), plt.yticks([])
	plt.show()
	def PSNREval(image1, image2, R):
	error = np.sum((image1.astype(float) - image2.astype(float)) ** 2)
	MSE = error / (float(image1.shape[0] * image1.shape[1]))
	PSNR = 10 * math.log10(255 * 255 / MSE)
	print("MSE: " + str(MSE))
	print("PSNR: " + str(PSNR))
	if __name__ == "__main__":
	main()

view raw ImageFFTFilter.py hosted with by GitHub

コードの解説

画像への高速フーリエ変換

	# FFT (time domain -> frequency domain)
	image_fq = np.fft.fft2(image_tm);

view raw ImageFFTFilter.py hosted with by GitHub

上記の箇所で画像に対してFFTを実施し、周波数成分へと変換します。

周波数成分のシフト

image_fq_shifted = np.fft.fftshift(image_fq)

view raw ImageFFTFilter.py hosted with by GitHub

デフォルトだと低周波成分が四隅に現れるので、fftshiftを用いて低周波成分が画像の中央に来るようにシフトします。

特定の周波数成分のマスク

	# Frequency mask creation
	CUTOFF_FREQ = 50;
	size = image_tm.shape
	mask = np.zeros(size)
	length_x = size[1]
	length_y = size[0]
	center_x = size[1]/2
	center_y = size[0]/2
	for x in range(0,length_x):
	for y in range(0,length_y):
	if abs(x - center_x) < CUTOFF_FREQ and abs(y - center_y) < CUTOFF_FREQ:
	mask[x,y]=1
	else:
	mask[x,y]=1e-10

view raw ImageFFTFilter.py hosted with by GitHub

上記の箇所で特定の周波数成分をマスクします。

「CUTOFF_FREQ = 50」が残す周波数成分の閾値を表しており、この値が大きいほど高周波成分まで残ります。

周波数成分のシフト（戻し）

image_fq = np.fft.fftshift(image_fq_shifted)

view raw ImageFFTFilter.py hosted with by GitHub

周波数成分をシフトさせていたのを元に戻します。

画像の高速フーリエ逆変換

	# IFFT (frequency domain -> time domain)
	image_tm_inverse = np.fft.ifft2(image_fq).real;

view raw ImageFFTFilter.py hosted with by GitHub

上記の箇所で周波数成分に対して逆変換を実施し、元の画像へと戻します。

実行結果

Frequency Filterの白の部分が通過させる周波数、黒の部分が遮断する周波数になります。

CUTOFF_FREQ = 10

わずかな低周波成分のみを通過させる場合は、かなり元の画像に劣化が加わっていることがわかります。

CUTOFF_FREQ = 50

一方、ある程度低周波成分を通過させた場合には、かなり元の画像の特徴を保っていることがわかります。

このことは、画像を周波数領域に変換し、低周波数成分のみを保存することで、画像の大部分の特徴を保持しながら、データサイズを削減できることを示します。

まとめ

本日は画像のフーリエ変換を行い、特定の周波数成分を遮断してみました。

低周波成分のみを残すことでノイズ除去をしたり、高周波成分のみを残すことでエッジを抽出したり、フーリエ変換はさまざまな用途に応用することができます。