Pythonで画像処理〜デコンボリュージョン〜
こんにちは。
今回は、画像の劣化(ボケ)を直すという作業をしてみました。
大まかにいうと、点拡がり関数(Point Spread Function, PSF)というボケの原因の関数を求めて(推定で決めて)、それをフーリエ変換したOptical Transfer Function(OTF)を用いるとボケが解消されるといった流れです。
前回のPythonシリーズはこちら↓
kabos.hatenablog.jp
画像の準備(作成とぼかし)
今回は、まず自分でジーメンススターという単純なパターンを作成しました。
ここはコードと画像だけ載せておきます。
#サイズを決める sh_up = (2005, 2005) ncycle = 60. #座標作成 limit_y = (sh_up[0]-1) // 2 limit_x = (sh_up[1]-1) // 2 y_up = np.arange(-limit_y, limit_y+1) x_up = np.arange(-limit_x, limit_x+1) xx_up, yy_up = np.meshgrid(x_up, y_up) #ジーメンススター作成 angles = np.arctan2(xx_up, yy_up) fields = angles % (2. * np.pi / ncycle) star_up = fields * np.pi / ncycle #ダウンサンプリング nbin = 5 shape_star = star_up.shape shape_star_binned = [shape_star[0] // nbin, shape_star[1] // nbin] star = np.reshape(star_up, (shape_star_binned[0], nbin, shape_star_binned[1], nbin)) star = np.mean(star, (1,3)) plt.figure(figsize=(15, 8)) plt.subplot(121) plt.title('Up-sampled (binary)') plt.imshow(star_up, interpolation='none', cmap='gray') plt.subplot(122) plt.title('Down-sampled (smooth)') plt.imshow(star, interpolation='none', cmap='gray')
この模様に名前がついている事を知りませんでした。右がダウンサンプリング後のものですが、確かに中心部分などのジャギジャギ感が薄れています。
次に、ガウシアンフィルタを用いて画像をぼかします。
sh = star.shape sigma = 3 #座標作成 lw = 4 * sigma gauss = np.zeros((2*lw+1, 2*lw+1)) x, y = np.arange(-lw, lw+1), np.arange(-lw, lw+1) xx, yy = np.meshgrid(x, y) #カーネルの作成 gauss = np.exp(-(xx**2 + yy**2) / (2* sigma**2)) #フィルタの適用 star_gauss = nd.convolve(star, gauss, mode='wrap')
このように、画像をぼかすことができました。
ちなみに、ガウシアンフィルタについては最初の記事にも記述があります。
デコンボリュージョン(復元)
それではボケを戻していきます。この操作をデコンボリュージョンというようです。
まずは、フーリエ変換を行うためにカーネルのサイズをジーメンススターのサイズに合わせます。
gauss_padded = np.zeros(sh) gauss_padded[sh[0]//2-lw:sh[0]//2+lw+1, sh[0]//2-lw:sh[0]//2+lw+1] = gauss gauss_padded = np.fft.ifftshift(gauss_padded)
そしてPSTをフーリエ変換することでOTFを求め、たたみ込みの逆演算を行うことでデコンボリュージョンできます。
今回はボケの原因のカーネルがgaussなので、先ほどフーリエ変換用に直したgauss_paddedをフーリエ変換するとOTFが求まります。
#OTFを求める gauss_otf = np.fft.fft2(gauss_padded) #デコンボリュージョン star_gauss_deconvolved = np.fft.ifft2(np.fft.fft2(star_gauss)/gauss_otf).real
このように、無事復元できているのがわかります。
補足
ちなみに、OTFの振幅をModulation Transfer Function (MTF)、位相をPhase Transfer Function (PTF)といいます。それらを求めたのが次の画像です。
gauss_mtf = np.abs(gauss_otf) gauss_ptf = np.angle(gauss_otf)
