デコンボリューション、ボケ復元、WienerフィルターそしてRichardson-Lucy法2 - 光学設計とその周辺、そしてたまに全く関係ないやつ

さて2回目ではもう１つの方法であるRichardson-Lucy法を紹介しよう. 1回目の記事の通りデコンボリューションとはどうにかしてPSFの逆関数PSF^-1を求めて,以下の数式のとおり元の画像Iを得るということだった.

$I(x)= PSF^{-1}＊O \tag{1} \label{1}$

この時登場するのは入力画像I(x), 出力画像（測定画像）O(x'), PSF(x,x')3つであるがこれらを確率の考え方を用いて解釈してみます. すると光子1つに対して以下の通りに考えられるはず.

I(x):光子が物体面上のxにある確率P(x)
O(x'):光子が像面上のx'にきた確率P(x')
PSF(x,x'):物体面で光子1つがxにあった時、それが光学系を通して像面x'に到着する確率

この中でPSF(x,x')というのは光子がもともとxにあるという前提の時の光がx'につくときであり、こういったことを条件付き確率と呼び確率論ではP(x'|x)という記号で表したりする.
さらに逆関数PSF^-1(x,x')はx'が起きた時に, xが起きる確率と考えることもでき, こちらはP(x|x')と表される. これまで登場した4つの確率の関係式が以下の俗にいうベイズの定理である.分母の変換は要はx'の前提であるxについて全ての場合で積分すればP(x')そのものが得られる, ということです.

$P(x|x')=\dfrac{P(x'|x) P(x)}{P(x')} =\dfrac{P(x'|x) P(x)}{\int P(x'|x)P(x) dx} \tag{2} \label{2}$

前回の記事で説明したようにノイズのせいで直接求めるのが難しい逆関数PSF^-1をベイズの定理を使ってうまく求めているのがこのRL法のエッセンスである.
もちろんベイズの定理は今回のような光学・画像に限った話ではなく、ある2つの事象の確率とお互いに条件となる場合の２つの確率に一般に成立する理論である. 今回の例のように原因と結果の関係を逆にすることができるのがベイズの定理が強力な理由です.

実際に最終的に求めたいのはP(x)である. $P(x)=\int P(x|x') P(x') dx$ となることを利用して(2)式から以下の通りP(x)が得られる.

$P(x)⇒\int P(x|x') P(x') dx' =\int \dfrac{P(x'|x) P(x) P(x')}{\int P(x'|x)P(x) dx} dx'=P(x) \int \dfrac{P(x'|x) P(x')}{\int P(x'|x)P(x) dx} dx' \tag{3} \label{3}$

元の記号I、O、PSFを使うと最終的には以下の式となる.

$I(x)⇒I(x) \int \dfrac{PSF(x-x') O(x') }{\int PSF(x'-x) I(x) dx} dx'=I(x) \{ [ \dfrac{O(x) }{ PSF(x)＊ I(x) }] ＊PSF(-x) \} \tag{4} \label{4}$

左端のPSF(x)が反転してPSF(-x)になることに注意. このままだと右辺にも求めたいI(x)が含まれているので反復式で表すとi+1番目のI(x)はi番目に計算したI(x)を使うと以下の通りになる.

$I_{i+1}(x)=I_{i}(x) \{ [ \dfrac{O(x)}{I_{i}(x)＊PSF(x) } ] ＊PSF(-x) \} \tag{5} \label{5}$

これを反復的に求めていくのがRL法である. RL法を説明する文献の中には唐突にベイズの定理が現れるものもあり混乱するかもしれないが, これまでの説明の通りベイズの定理そのものと画像に関する諸量がどのように確率として解釈できるかがわかればそこまで難しい話ではないと思う. 前回紹介したデコンボリューションの方法と比べると, RL法は周波数領域に変換していないため, ややこしい発散の問題は発生しない.

せっかくなのでpythonを使って実装してみよう. 既にskimageでRL法のライブラリが用意されているしサンプルコードもある.
scikit-image.org

同じことをしても芸がないので, そのサンプルコードを参考に自分で実装してみました.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import convolve2d as conv2
from skimage import data, restoration

image_raw = data.camera()
image_original=np.array(image_raw/255,dtype=np.float) #1に正規化

# psfを定義して画像をぼかす
psf = np.ones((5, 5))/25
image_blurred = conv2(image_original, psf, 'same')

# ノイズを加える
image_noisy = image_blurred.copy()
image_noisy += (np.random.poisson(lam=25, size=image_raw.shape) - 10) / 200.

# ノイズありのボケた画像をRL法で復元
deconvolved_RL = restoration.richardson_lucy(image_noisy, psf, iterations=10)

fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(5, 5))
plt.gray()

for a in (ax[0,0], ax[0,1], ax[1,0],ax[1,1]):
       a.axis('off')

ax[0,0].imshow(image_original)
ax[0,0].set_title('image without blur')

ax[0,1].imshow(image_blurred)
ax[0,1].set_title('blurred image')

ax[1,0].imshow(image_noisy)
ax[1,0].set_title('blurred image with noise')

ax[1,1].imshow(deconvolved_RL, vmin=image_noisy.min(), vmax=image_noisy.max())
ax[1,1].set_title('RL image')

fig.subplots_adjust(wspace=0.02, hspace=0.2,top=0.9, bottom=0.05, left=0, right=1)
plt.show()

結果は以下の通り. RL法で復元した画像(右下)は元のボケ, ノイズ無しの左上の画像に近づいている.

ちなみに元のソースを見ると以下の通り実装されている. psf_mirror = np.flip(psf)が式(5)の右端の座標をマイナスにしているの対応しており, 非常にシンプルであることがソースを見ても分かる.

#https://github.com/scikit-image/scikit-image/blob/master/skimage/restoration/deconvolution.pyより一部改変

def richardson_lucy(image, psf, iterations=50, clip=True, filter_epsilon=None):
    float_type = np.promote_types(image.dtype, np.float32)
    image = image.astype(float_type, copy=False)
    psf = psf.astype(float_type, copy=False)
    im_deconv = np.full(image.shape, 0.5, dtype=float_type)
    psf_mirror = np.flip(psf)

    for _ in range(iterations):
        conv = convolve(im_deconv, psf, mode='same')
        if filter_epsilon:
            relative_blur = np.where(conv < filter_epsilon, 0, image / conv)
        else:
            relative_blur = image / conv
        im_deconv *= convolve(relative_blur, psf_mirror, mode='same')

    if clip:
        im_deconv[im_deconv > 1] = 1
        im_deconv[im_deconv < -1] = -1

    return im_deconv