チーム20B3

課題名

空間・周波数フィルタリング(Wiener filter, LPF, Gaussian filter)

研究者名

Motoki Takematsu
Takahiro Kawano

概要

テストのためのノイズ画像生成プログラムと、
Wiener filter、LPF、Gaussian filterを作成した。
実装はPython3.6.5、実行はUbuntu上で行った。

ソースコード

import numpy as np
from numpy.fft import fftshift, fft2, ifft2
from scipy.signal import gaussian, convolve2d
import matplotlib.pyplot as plt

#ブラー画像の生成
def blur(img, kernel_size = 3):
  h = np.eye(kernel_size) / kernel_size
   
   dummy = convolve2d(img, h, mode = 'same')
   return np.uint8(dummy.real)

#ガウスカーネルの生成
def gaussian_kernel(kernel_size = 3):
  h = gaussian(kernel_size, kernel_size / 3).reshape(kernel_size, 1)
   h = np.dot(h, h.transpose())
   h /= np.sum(h)
   return h

#ガウシアンノイズの付与
def add_gaussian_noise(img, sigma):
  gauss = np.random.normal(0, sigma, np.shape(img))
   noisy_img = img + gauss
   noisy_img[noisy_img < 0] = 0
   noisy_img[noisy_img > 255] = 255

  return noisy_img

#画像をグレースケール化
def rgb2gray(rgb):
  return np.dot(rgb[...,:3], [0.2989, 0.5870, 0.1140])

#Lowpassフィルタ
def low_pass_filter(img, r=0.3):
  src = np.copy(img)
   src = fft2(src)

  #画像の中心と、採用する範囲を決定
   cx, cy = int(img.shape[1]/2), int(img.shape[0]/2)
   rx, ry = int(r*cx), int(r*cy)

  fsrc = fftshift(src)
   fdummy = np.zeros(src.shape, dtype = complex)
   fdummy[cy-ry:cy+ry, cx-rx:cx+rx] = fsrc[cy-ry:cy+ry, cx-rx:cx+rx]
   fdummy2 = fftshift(fdummy)
   dummy = ifft2(fdummy2)

  return np.uint8(dummy.real)

#gaussianフィルタ
def gaussian_filter(img, kernel):
  return convolve2d(img, kernel, mode = 'same')

#wienerフィルタ及びinverseフィルタ
def wiener_filter(img, kernel):
  kernel /= np.sum(kernel)
   dummy = np.copy(img)
   dummy = fft2(dummy)
   dummyfft = 20*np.log(np.abs(dummy))
   kernel = fft2(kernel, s = img.shape)
   
   #wienerフィルタ
   wiener = np.conj(kernel) / (np.abs(kernel) ** 2+0.01)
   wienerfft = 20*np.log(np.abs(wiener))
   dummy1 = dummy * wiener
   dummy1 = np.abs(ifft2(dummy1))
   
   #inverseフィルタ
   inv = np.conj(kernel) / (np.abs(kernel) ** 2)
   invfft = 20*np.log(np.abs(inv))
   dummy2 = dummy * inv
   dummy2 = np.abs(ifft2(dummy2))
   
   plt.subplot(131),plt.imshow(dummyfft, cmap = 'gray')
   plt.title('blur fft'), plt.xticks([]), plt.yticks([])
   plt.subplot(132),plt.imshow(wienerfft, cmap = 'gray')
   plt.title('wiener fft'), plt.xticks([]), plt.yticks([])
   plt.subplot(133),plt.imshow(invfft, cmap = 'gray')
   plt.title('inverse fft'), plt.xticks([]), plt.yticks([])
   
   plt.show()
   return dummy1, dummy2

#画像間の類似度の尺度であるPSNRの計測
def psnr(img_1, img_2, data_range=255):
  mse = np.mean((img_1.astype(float) - img_2.astype(float)) ** 2)
   return 10 * np.log10((data_range ** 2) / mse)

if __name__ == '__main__':
  #グレースケールで画像の読み込み 
   img = rgb2gray(plt.imread('lena.jpg'))
   
   #ブラー画像の生成
   blurred_img = blur(img, kernel_size = 9)
   
   #ノイズ画像の生成
   noisy_img = add_gaussian_noise(img, sigma = 30)

  #Lowpassフィルタ処理
   low_pass = low_pass_filter(noisy_img, r=0.5)
   
   #wienerフィルタ処理及びinverseフィルタ処理
   f_kernel = np.eye(9) / 9
   wiener, inv = wiener_filter(blurred_img, f_kernel)

  #gaussianフィルタ処理
   kernel = gaussian_kernel(9)    
   gauss = gaussian_filter(noisy_img, kernel)
       
   # Display results
   display = [img, blurred_img, inv, wiener]
   label = ['Original', 'Blur', 'Inverse Filter', 'Wiener Filter']
   
   fig = plt.figure(figsize=(12, 10))
       
   for i in range(len(display)):
       fig.add_subplot(2, 2, i+1)
       plt.imshow(display[i], cmap = 'gray')
       plt.title(label[i])
   plt.show()

  display = [img, noisy_img, low_pass, gauss]
   label = ['Original', 'noise', 'low pass', 'gaussian']
   psnrlist = ['', psnr(img, noisy_img), psnr(img, low_pass), psnr(img, gauss)]
   
   fig = plt.figure(figsize=(12, 10))
       
   for i in range(len(display)):
       fig.add_subplot(2, 2, i+1)
       plt.imshow(display[i], cmap = 'gray')
       plt.title(label[i])
       plt.xlabel(psnrlist[i])
       
   plt.show()

テスト画像の生成

9×9のモーションフィルタを畳み込むことで、ブラー画像を生成した。
Figure2_2re.jpg
平均0、標準偏差30のガウシアンノイズを加えることで、ノイズを含む画像を生成した。
Figure3_2.jpg

実行結果

Figure_fft.png
左がブラー画像のFFT結果、真ん中がwienerフィルタのFFT結果、右がwienerフィルタの簡易版であるinverseフィルタのFFT結果である。
各周波数スペクトルの積をとることで、周波数フィルタリングが可能である。

Figure_kekka1.jpg
左下はinverseフィルタ処理、右下はwienerフィルタ処理を行った結果である。
wienerフィルタで用いる劣化関数は9×9のモーションフィルタを用いている。wienerフィルタリングによって、元画像の復元がかなりできていることがわかる。

Figure_kekka2.jpg
左下はLPF処理、右下はガウシアンフィルタ処理を行った結果である。
LPFに関しては、ノイズ画像を高速フーリエ変換した後、ゼロ周波数シフトをかけて低周波数成分を中心に集め、FFT画像の中心から半径rの部分だけを採用することで実装した。ガウシアンフィルタは、9×9のガウスカーネルとの畳み込み演算で実装した。
どちらもノイズ除去としての役割があるが、見た目ではどちらが優れているかわかりづらいため、画像の類似度の尺度として用いられるPSNRを実装し、元画像との類似度を計測している。



考察

Inverseフィルタでは予期せぬ結果が得られてしまったが、これは元々Inverseフィルタがノイズに敏感であり、畳み込みの際に外周部の0を含めた計算を行ってしまったためであると考えられる。
またWienerフィルタではノイズ対信号比を考慮することでノイズの影響を軽減することができたが、リンギングと呼ばれる波型の模様が生じた。
このリンギングは画像復元において深刻な課題であり、リンギングを除去できるようなフィルタリングが盛んに研究されている。
また、劣化関数も劣化画像を生成したものと同様のものを用いているため、うまく復元できるのは当然かもしれない。実際にはこの劣化関数やノイズ対信号比の値は未知であるため、画像の復元は容易な課題ではない。
LPFとGaussianフィルタはほぼ同レベルで元画像を再現することができたが、PSNRの値を見るとGaussianフィルタの方が上手くいったと考えられる。

今後の展望として、ブラー画像を生成するときに、convolve2dのモードを変更することで逆フィルタがうまく適用できる場合もあるようなので、そちらも試してみたい。また、今回は試していないが、ブラーとノイズの両方をかけた劣化画像がwienerフィルタによってうまく復元できるか試してみたい。ノイズ対信号比の推定アルゴリズムもあるようなので実装してみるとよいかもしれない。

参考文献


  • 最終更新:2020-06-30 17:14:11

このWIKIを編集するにはパスワード入力が必要です

認証パスワード