フーリエ変換に基づき画像の空間周波数を確認する（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 image data
        
        
          
              filename = "data/noisy_image.png"
        
        
          
              image = cv2.imread(filename, cv2.IMREAD_GRAYSCALE)
        
        
          
              
        
        
          
              if image is None:
        
        
          
                  print("Cannot find image : " + filename)
        
        
          
                  sys.exit()
        
        
          
          

        
        
          
              image_tm = np.asarray(image)
        
        
          
              
        
        
          
              # FFT (time domain -> frequency domain)
        
        
          
              image_fq = np.fft.fft2(image_tm);
        
        
          
              # IFFT (frequency domain -> time domain)
        
        
          
              image_tm_inverse = np.fft.ifft2(image_fq).real;
        
        
          
              
        
        
          
              # Component replacement (Low-frequency components -> center)
        
        
          
              image_fq_shifted = np.fft.fftshift(image_fq)
        
        
          
              # Calculating magnitude spectrum
        
        
          
              magnitude_spectrum = 20 * np.log(np.abs(image_fq_shifted));
        
        
          
          

        
        
          
              # Evaluation with PSNR
        
        
          
              PSNREval(image_tm, image_tm_inverse, 255)
        
        
          
          

        
        
          
              plt.subplot(131),plt.imshow(image_tm, cmap = 'gray')
        
        
          
              plt.title('Image (input)'), plt.xticks([]), plt.yticks([])
        
        
          
              plt.subplot(132),plt.imshow(magnitude_spectrum, cmap = 'gray')
        
        
          
              plt.title('Magnitude spectrum'), plt.xticks([]), plt.yticks([])
        
        
          
              plt.subplot(133),plt.imshow(image_tm_inverse, cmap = 'gray')
        
        
          
              plt.title('Image (inverse)'), 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
        
          ImageFFT.py
        
        hosted with  by GitHub

コードの解説

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


    
      
        
  
    
    

        


  


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



    

  


      
      
        view raw
        
          ImageFFT.py
        
        hosted with  by GitHub

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

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


    
      
        
  
    
    

        


  


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



    

  


      
      
        view raw
        
          ImageFFT.py
        
        hosted with  by GitHub

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

周波数成分の大きさの確認


    
      
        
  
    
    

        


  


  
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
          
              # Component replacement (Low-frequency components -> center)
        
        
          
              image_fq_shifted = np.fft.fftshift(image_fq)
        
        
          
              # Calculating magnitude spectrum
        
        
          
              magnitude_spectrum = 20 * np.log(np.abs(image_fq_shifted));
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
  



    

  


      
      
        view raw
        
          ImageFFT.py
        
        hosted with  by GitHub

上記は確認用です。周波数成分の大きさを確認します。ここで、何もしないと低周波成分が四隅に現れるので、fftshiftを用いて低周波成分の大きさが画像の中央に来るようにシフトします。

実行結果

以下のような結果が得られます。真ん中の画像で、白い部分が成分が多いことを示します。中央に近づくほど低周波成分なので、低周波成分が多めであることがわかります。

また、今回は周波数成分に変化を与えずに元に戻しているため、「Image(input)」と「Image(inverse)」は一致します。

本当に一致しているかを確認するために、コンソールにMSEとPSNRを出力できるようにしています。

実行したところ「MSE: 9.393230390470507e-28」と、ほんのわずかに差異があるようですが、浮動小数点の誤差などでしょうか。ただ、この数値であれば一致していると言って良いかと思います。

MSEとPSNRは、2つの画像がどれくらい一致しているかを示す評価指標です。以下の記事で紹介しているので、興味のある方はご確認ください。

まとめ

本日は画像のフーリエ変換を紹介しました。次回はフーリエ変換で得た周波数成分に何か変化を加えて、周波数成分を使った画像処理を行ってみましょう。

以下の記事で特定の周波数を除去するプログラムも紹介していますので、参考にしてみてください。

	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 image data
	filename = "data/noisy_image.png"
	image = cv2.imread(filename, cv2.IMREAD_GRAYSCALE)

	if image is None:
	print("Cannot find image : " + filename)
	sys.exit()

	image_tm = np.asarray(image)

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

	# Component replacement (Low-frequency components -> center)
	image_fq_shifted = np.fft.fftshift(image_fq)
	# Calculating magnitude spectrum
	magnitude_spectrum = 20 * np.log(np.abs(image_fq_shifted));

	# Evaluation with PSNR
	PSNREval(image_tm, image_tm_inverse, 255)

	plt.subplot(131),plt.imshow(image_tm, cmap = 'gray')
	plt.title('Image (input)'), plt.xticks([]), plt.yticks([])
	plt.subplot(132),plt.imshow(magnitude_spectrum, cmap = 'gray')
	plt.title('Magnitude spectrum'), plt.xticks([]), plt.yticks([])
	plt.subplot(133),plt.imshow(image_tm_inverse, cmap = 'gray')
	plt.title('Image (inverse)'), 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()