はじめに
Richardson-Lucy deconvolution（以下、RL法）は、画像等を鮮明にする手法である。一般的な撮像観測では、観測装置の回折等により光源が拡がって撮像される。この拡がり具合は、点拡がり関数(Point Spread Function, 以下、PSF)で表される。RL法は既知のPSFと撮像画像から、真の画像を推定する手法である。本記事ではRL法の理論式を導出する。
RL法の理論式
RL法は、既知のPSFと撮像画像を用いて真の画像を推定する手法である。真の画像の分布、既知のPSF、撮像画像の関係をベイズの定理を用いて関係づける。そして、ベイズ推定を繰り返し計算することにより真の画像を推定する。理論式は
W_{i,r+1} = W_{i,r}\sum_k \frac{P_{i,k}H_k}{\sum_j P_{j,k}W_{j,r}}\
\quad r = 0,\ 1,\ 2,\ \dots \tag{1}
と表される。$i,j,k$はピクセル座標とする。$W$は真の画像、$H$は撮像画像、$r$は反復回数とする。$P_{j,k}$はPSFで、$W_j$が$H_k$で観測される確率を$P(H_k|W_j)=P_{j,k}$とする。
RL法の理論式の導出
式(1)を示す。$W_i$を真の画像の$i$番目のピクセル値、$H_k$を撮像画像の$k$番目のピクセル値とする。以下、ピクセル値では分かりにくいため、光子として考える。
任意の1つ光子が$H_k$であるとき、$W_i$からのものである確率は、ベイズの定理で
P(W_i|H_k) =
\frac{P(H_k|W_i)P(W_i)}{\sum_{j} P(H_k|W_j)P(W_j)} \tag{2}
と表される。まず、式(2)の左辺と右辺の意味を考えることで、式(2)全体の意味を理解する。
- 左辺は、$H_k$に入った全ての光子のうち、$W_i$からきた光子の確率を表す。
- 右辺の分子は、任意の1つの光子が$W_i$でもあり$H_k$でもある確率を表す。
- 右辺の分母は、任意の1つの光子が$H_k$である確率を表す。（シグマ和を取るところ以外は、右辺の分子同様に任意の1つの光子が$W_j$でもあり$H_k$でもある確率を表す。それを全てのピクセル座標で和を取るため、真の画像$W$の任意のピクセル座標のものからでよい。つまり、任意の1つの光子が$H_k$である確率を表す。）
- 右辺の分子分母から右辺は、$H_k$に入った全ての光子のうち、$W_i$からきた光子の確率を表す。
よって、式(2)の左辺と右辺の同じ意味である。このように、ベイズの定理によって求めたい関係性を既知の情報を含んだ別の関係性から捉えることができる。また、式(2)の$P(H_k|W_i)$は尤度と呼ばれ、既知のPSFで表される。
ここで、$W_i$がPSF程度に散らばって$H_k$で得られるので、全ての$W_i\cap H_k$を計算することで$W_i$を求めると、
P(W_i) =
\sum_{k} P(W_i\cap H_k)=\sum_{k} P(W_i|H_k)P(H_k) \tag{3}
となる。（条件付き確率より、$P(W_i|H_k)=P(W_i\cap H_k)/P(H_k)$を用いた。）
式(2)を式(3)に代入すると、
P(W_i) =
\sum_{k} \frac{P(H_k|W_i)P(W_i)P(H_k)}{\sum_{j} P(H_k|W_j)P(W_j)} \tag{4}
となる。
$P(W_i)$は本来知り得ない値なので、適当な推定値（撮像画像、全て均等な値の画像等）を初期値として入れることで逐次的に求めることができ、
P_{r+1}(W_i) =
P_r(W_i)\sum_{k} \frac{P(H_k|W_i)P(H_k)}{\sum_{j} P(H_k|W_j)P_r(W_j)} \quad
r = 0,\ 1,\ 2,\ \dots \tag{5}
となる。
PSFが規格化されていて全てのピクセル座標で同一のものを使用するとき、$W_i$のピクセルの座標を中心として$H_k$の確率は、
P(H_k|W_i) = P_{i,k} \tag{6}
となる。また、真の画像と撮像画像で光子の総数は不変なので、
\begin{align}
P(W_i) &= \frac{W_i}{\sum_l W_l} \\
P(H_k) &= \frac{H_k}{\sum_l H_l} = \frac{H_k}{\sum_l W_l} \tag{7}
\end{align}
となる。
よって、式(6), (7)を式(5)に代入して
\begin{align}
\frac{W_{i,r+1}}{\sum_l W_l} &=
\frac{W_{i,r}}{\sum_l W_l}
\sum_{k} \frac{P_{i,k} \frac{H_k}{\sum_l W_l}}
{\sum_j P_{j,k}\frac{W_{j,r}}{\sum_l W_l}} \\
\iff W_{i,r+1} &= W_{i,r}\sum_k \frac{P_{i,k}H_k}{\sum_j P_{j,k}W_{j,r}}
\quad r = 0,\ 1,\ 2,\ \dots \tag{9}
\end{align}
となる。よって、RL法の理論式(1)は示された。
参考文献
- Richardson, William Hadley. 1972, JOSA, 62, 55–59
- Lucy, L. B. 1974, Astronomical Journal, 79, 745–754