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

2回分の記事を使って光学とは少し離れてちょっと画像処理よりのことを取り上げてみます. 俗にいうデコンボリューションとかボケ復元とか言われるトピックです. 特に2回目の記事はRichardson-Lucy法(RL法)をメインに取り上げます. よくあるRL法の説明ではやや唐突気味に条件付き確率, ベイズ推定などの頭が痛くなりそうなキーワードが登場し, 自分もそうでしたが混乱するかもしれないです. ただ実際にRL法の導出でベイズ推定の理論が使われるのはたった1つの式のみで, しかもそれ自体は初等的な確率論に少し毛が生えた程度です. よって決して難しい話ではなくむしろRL法を考えた人は(その名の通りRichardsonとLucy)はうまいこと考えたな, と本記事を読んだ人に感じてもらえば幸いです.

まず第１回目の記事ではデコンボリューションの概要とWeinerフィルターと呼ばれる方法を取り上げてみよう.

デコンボリューションの概要
単純なデコンボリューション
Wienerフィルター

デコンボリューションの概要

レンズのような光学系を通した像の分布は, 元の画像の分布 $I(x)$ と光学系を通した像分布 $O(x)$ , そしてその光学系のボケ具合を表した $PSF(x)$ を使って以下の式で表される. ＊は畳み込み計算の記号.

$\displaystyle O(x)=\int PSF(x-x')I(x') dx' = PSF＊I \tag{1} \label{1}$

もし光学系がボケのない理想的なレンズ系なら $PSF$ は $\delta$ 関数となり、像がそのまま元の画像Iとして取得できる. ボケ復元とは直接の測定値である $O(x)$ と既知である $PSF(x)$ の情報からボケのない元の画像 $I(x)$ を得ることだ. そのアプローチとして $PSF$ の逆関数 $PSF^{-1}$ (MTFのこと)を何とかして求めて $O(x)$ に計算し, $I(x)$ を求めることをデコンボリューションと呼ぶ. 書くほどではないが式で表すと以下の通り.

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

この逆関数 $PSF^{-1}$ を求めるのに古今東西いろんな手法が提案されてきた.

単純なデコンボリューション

もしフーリエ解析を学んだことがあれば「えっ？こんなの式(1)のフーリエ変換をすれば右式は単に２つの値の積になるから, $PSF$ のフーリエ変換で両辺を割って, 再度逆フーリエ変換すればいいでしょ？」となるかもしれない. つまり以下の式のとおり表すことです.

$\mathcal{F} [O]=\mathcal{F} [PSF] \mathcal{F}[I] \tag{3} \label{3}$

$\displaystyle I= \mathcal{F}^{-1} [ \dfrac {\mathcal{F}[O]}{\mathcal{F}[PSF] } ] \tag{4} \label{4}$

ここで $1/{\mathcal{F}[PSF] }$ の部分を逆フィルタと呼ぶ.
これはもちろん間違いではないですが, 実際にこの計算を実行してもうまくいかないことが多い.
というのは実際の測定値にはノイズが含まれており, (4)の計算過程でこのノイズが強調され良い画像が得られないからだ. 具体的には $\mathcal{F} [PSF]$ が0に近いような周波数のところではそこで式(4)が発散し, その空間周波数のノイズが増幅されてしまうなんてことが起きる.

Wienerフィルター

とにかくノイズの増幅を抑えることがデコンボリューションでは重要となる.
改めて式(1)を以下の通りに書き直してみる.

$\displaystyle \displaystyle O(x)=\int PSF(x-x')I(x') dx' = PSF＊I +n(x) \tag{5} \label{5}$

ここで $n(x)$ とはノイズ成分を表す. (3)式の逆フィルタの部分を $\mathcal{F}[PSF]$ を $H$ とおいて以下の式で置き換えてみよう.

$\displaystyle \displaystyle \dfrac{1}{H}⇒\dfrac{1}{H } \dfrac{|H|^{2}}{|H|^{2} +| \mathcal{F}[ n]|^{2}/|\mathcal{F}[ I]|^{2} } \tag{6} \label{6}$

このように元画像で正規化して考えておくと発散を抑えることができる.

といっても肝心のその元画像 $I(x)$ とノイズ $n(x)$ がわからないので, 実際に実装するときは(6)式の該当の部分を適当なパラメーター $\Gamma$ で置き換えて以下のようにする.

$\displaystyle \dfrac{1}{H}⇒\dfrac{1}{H } \dfrac{|H|^{2}}{|H|^{2} + \Gamma } \tag{7} \label{7}$

対処療法的なところがあるが, このようにして逆関数を表すことをWienerデコンボリューションと呼び, 式(7)のフィルターをWeinerフィルターと呼ぶことが多い.

以上で１回目は終わり. 2回目は別のRichardson-Lucy法を扱います.
eikonal.hatenablog.jp