はじめに
本記事では、天文学などで非常に有名な古典的手法のRichardson-Lucy deconvolution (RL法) について、はじめから丁寧に説明したいと思う。私自身、RL法とベイズ推定や最尤法などの関係性をあまり理解できていないが、これまでの研究（学部4年~修士1年の約2年間）で得られた範囲で、できるだけ丁寧にまとめたつもりである。本記事により、RL法を使ってみたい初学者や、RL法を知っているけどもっと理解を深めたい方の一助になればと思う。
ブレとは
ブレ（ブラー、ぼけ）とは、画像が手ブレ、カメラの動き、被写体の動きなどにより不鮮明に見える現象である。以下が典型的なブレ画像の一例である。
左図の鮮明な画像に比べて、右図の画像全体がぼんやりとした画像がブレ画像と呼ばれる。このブレ画像の数学的に理解するとき、点源を撮影したときにどのように広がって観測されるかを表した、点拡がり関数（Point Spread Function, 以下、PSF）と呼ばれるものが使われる。上図の場合のPSFを入れると、以下のようになる。
このように、PSFで点がどのように広がるかがわかると、ある現実の複雑な画像においても観測時の不鮮明な画像が得られるプロセスについて知ることができる。その観測のプロセスを簡単に書くと、どんなに複雑な画像でも鮮明な画像の各ピクセルが点源だとして、それぞれがPSFのように広がって観測され、全てが重なり合ってブレとして観測される。これは、数学的には畳み込みで書くことができ、
H(x)=\int W(\xi)P(x|\xi)d\xi \tag{1}
と表される。ただし、この式は画像内で一様でかつノイズがない場合のものである。ここで、$W(\xi)$は知りたい元の鮮明な真の画像の$\xi$番目のピクセル値、$H(x)$は観測画像の$x$番目のピクセル値である。$P(x|\xi)$は、真の画像$W$の$\xi$番目の値が、観測画像$H$の$x$番目で観測される確率を表し、いわゆるPSFである。（本来は、式(1)で積分の記号を使っているので、$W(\xi)$は、真の画像の無限に小さな区間区間$[\xi,\xi+d\xi]$を意味することになるが、本記事では本質的ではないため、特に意識して書いていない。本質的な話を現実的な画像で置き換えると、「無限に小さいピクセルでできた無限に大きな画像を使っている」といったイメージである。）
ここで、PSFは確率分布のため、観測画像のあるピクセル$x$に入ってくる確率の和と、真の画像のある$\xi$から放射された確率の和は、
\int P(x|\xi)d\xi=\int P(x|\xi)dx=1 \tag{2}
を満たし、したがって、式(1)を積分すると
\begin{aligned}
\int H(x)dx&=\int \int W(\xi)P(x|\xi)d\xi dx\\
&=\int W(\xi)\left(\int P(x|\xi)dx\right) d\xi\\
&=\int W(\xi) d\xi\\
\end{aligned}
\tag{3}
となる。つまり、真の画像$W$と観測画像$H$の値の合計は保存されるというごく当たり前の事が確認できる。(式(3)の理解として例えば、カメラで撮影された画像が、手ブレした場合でもその明るさの総量は手ブレしなかった場合と同じであるイメージを持つとよいと思う。ただし、サチってしまうピクセルやブレによって撮像時の感度に違いでノイズが入り、総量が合わない場合はその限りではない。そういう意味では、現実の世界ではPSFが必ずしも1ではないが、ここは理想的な世界のためこのまま進める。）
真の画像全ての座標に対して、PSFの確率をかけて拡散させ足し上げるプロセスが観測であるが、数学的には畳み込み計算を行なっているため、畳み込み演算子$\otimes$を使って、
H(x)\stackrel{\mathrm{def}}{=}(W \otimes P)(x)\stackrel{\mathrm{def}}{=}\int W(\xi)P(x-\xi)d\xi \tag{4}\\
と定義されることもあり、本記事もこの式を一部使用する。
ここまでの話を踏まえ、観測のプロセスをイメージで表すと以下のようになる。
真の鮮明な画像の各ピクセルがPSFで畳み込まれて観測されることがわかるイメージ図である。
イメージデコンボリューションとは
イメージデコンボリョーションとは、ブレの原因となるPSFなどを手掛かりに、元の鮮明な画像を復元するための手法である。
イメージデコンボリューションの種類
PSFが既知、未知の場合でそれぞれ、Non-Blind deconvolution、Blind deconvolutionと呼ばれる。RL法は、何かしらの方法で得たPSFを使って真の画像を推定するため、Non-Blind deconvolutionに属する。ここで、イメージを使って説明する。
- 左図は、観測のイメージを表しており、真の画像がPSFで畳み込まれて、観測される様子を図で示したものである。
- 中央図は、Non-Blind deconvolutionのイメージで、既知のPSFと観測画像から逆畳み込みで真の画像を推定するイメージである。
- 右図は、Blind deconovlutionのイメージで、観測画像のみから真の画像を推定する手法である。ただし、真の画像を得ることが目的のため、必ずしも未知のPSFを推定しないで求める手法もある。
RL法とは
RL法とは、ベイズ推定によって反復的に真の鮮明な画像を推定する手法であり、
W_{r+1}(k)=W_r(k)\left(\frac {H}{W_r \otimes P}\otimes P^*\right)(k)\ \quad r = 0,\ 1,\ 2,\ \dots \tag{5}
と表せる。ここで、$\otimes$は畳み込み演算子、$k$は画像の各画素に対応した座標、$r$は反復回数、$W$は元の鮮明な真の画像、$P$はPSF、$P^*$はPSFを各軸に対して反転した配列(特に2次元画像の場合は$P$を180度回転した画像)、$H$は撮像された画像である。$W_{0}$に適当な初期値（一般的に全ての画素に同じ値を入れた画像が使われることが多い)を入れてこの手法を始める。反復処理により収束した真の画像は、最尤解に近づくことが知られている。
まず、RL法を適用した例を以下に示す。
以下は、RL法の反復回数1回から80回までの結果である。
この例では、PSFで畳み込み観測画像を用意した。そして、畳み込んだPSFと同じもので、真の画像の初期値$W_0$の全ての画素を同じ値を入れたものを使い、RL法を行った。このように、数回の反復で急激に真の画像に近づき、収束していく様子がわかる。
このRL法の式の導出過程について以下で詳しく説明する。
観測とベイズの定理の関係
私たちが知りたい真の画像$W$を求めるためには、観測画像$H$の$x$の画素での値が、真の画像$W$の$\xi$の画素から来た確率$Q(\xi|x)$を計算する必要がある。この$Q(\xi|x)$は、観測のPSFの逆畳み込みの計算を意味するため、感覚的に理解することが難しいと思う。というのも、あるピクセルで観測される値は、真の画像の周辺の画素がPSFで畳み込まれて入ってくるため、観測画像から単純にどこのピクセルからどの程度やって来たのかを推測するのは真の画像の分布を知っていないとわからないと予想できるからである。イメージでも補足すると以下のようになる。
逆畳み込み$Q(\xi|x)$の計算が複雑であることがイメージ図からもわかると思う。この$Q(\xi|x)$は、数学的には、ベイズの定理の条件付き確率で記述することができ、
Q(\xi|x) = \frac{W(\xi)P(x|\xi)}{\int W(\xi)P(x|\xi)d\xi} \tag{6}
と表される。この関係式により、左辺の難しいPSFの逆畳み込み計算を右辺の既知のPSFを使った関係$P(x|\xi)$（尤度と呼ばれる）で、表すことができる。
ここで、右辺各項の意味を説明する。
- 右辺の分母 ($\int W(\xi)P(x|\xi)d\xi$): 真の画像$W$のPSFで畳み込み（観測のプロセス）の計算をしており、すなわち、観測画像$H(x)$を意味する。
- 右辺の分子 ($W(\xi)P(x|\xi)$): 真の画像$W$の$\xi$番目の画素値が、観測画像$H$の$x$番目の画素で観測される時の値を計算することを意味する。
畳み込み$P(x|\xi)$と、逆畳み込み$Q(\xi|x)$の間には、以下が成り立つ。
\begin{equation*}
\begin{cases}
H(x)=\int W(\xi)P(x|\xi)d\xi \ (畳み込み)\\
W(\xi)=\int H(x)Q(\xi|x)dx \ (逆畳み込み) \tag{7}
\end{cases}
\end{equation*}
式(7)の畳み込み計算は自明だが、逆畳み込み計算については、ベイズの定理で得られた数式を計算して確認する。
\begin{aligned}
\int H(x)Q(\xi|x)dx &= \int H(x)\frac{W(\xi)P(x|\xi)}{\int W(\xi)P(x|\xi)d\xi} dx\\
&= \int H(x)\frac{W(\xi)P(x|\xi)}{H(x)} dx\\
&= \int W(\xi)P(x|\xi) dx\\
&= W(\xi)\int P(x|\xi) dx\\
&= W(\xi) \ \left(\because \int P(x|\xi)dx=1\right)\\
\end{aligned}
\tag{8}
となり、式(7)の逆畳み込みの関係式は成り立つ。
ベイズ推定とRL法の関係
ベイズ推定とは、初期の事前分布(真の画像)を設定し、データを観測して事後分布(更新された真の画像)を計算し、次の事前分布として使うというプロセスを複数回繰り返すことで、真のパラメータを推定していく手法である。
ここで、ベイズ推定とRL法の関係について説明する。真の画像は本来知り得ないので、適当な主観的な$W_r$をベイズの定理の式(6)の$W$に代入すると、逆畳み込みの確率$Q_r$は、
Q_r(\xi|x) = \frac{W_r(\xi)P(x|\xi)}{\int W_r(\xi)P(x|\xi)d\xi} \tag{9}
と表される。事前分布$W_r$で求めた$Q_r$の関係式を式(7)の2つ目の逆畳み込みの式に代入すると、事後確率$W_{r+1}$は、
\begin{aligned}
W_{r+1}(\xi)&=\int H(x)Q_r(\xi|x)dx\\
&=\int H(x) \frac{W_r(\xi)P(x|\xi)}{\int W_r(\xi)P(x|\xi)d\xi}dx\\
&=W_r(\xi)\int H(x) \frac{P(x|\xi)}{\int W_r(\xi)P(x|\xi)d\xi}dx\\
&=W_r(\xi)\int \frac{H(x)}{H_r(x)}P(x|\xi)dx \ \left(H_r\colon =\int W_r(\xi)P(x|\xi)d\xi\right)
\end{aligned}
\tag{10}
となる。ここで、$H_r$は、$W_r$をPSFで畳み込んで得られた擬似的な観測画像を表す。式(10)の初期値$W_0$を決めてあげると、反復的に真の画像を推定する。反復のはじまりを$r=0$として書き下すと
W_{r+1}(\xi)=W_r(\xi)\int \frac{H(x)}{H_r(x)}P(x|\xi)dx\ \quad r = 0,\ 1,\ 2,\ \dots \tag{11}
となる。この式(11)はまさにRL法の式であり、式(5)と同等の計算であることも確認できる。
RL法の非負性
$W_0>0$のとき、式(11)の全ての項($H, H_0, P(x|\xi)$)も0より大きいので、$W_{r+1}>0$が成り立つ。このように、RL法で推定される真の画像には非負性の性質がある。
RL法の保存性
RL法は、真の画像と観測画像の積分値が保存する性質がある。この特徴について、式を追って説明する。
式(11)の両辺を真の画像全体で積分すると、
\begin{aligned}
\int W_{r+1}(\xi)d\xi&=\int \left(W_r(\xi)\int \frac{H(x)}{H_r(x)}P(x|\xi)dx \right) d\xi\\
&=\int \frac{H(x)}{H_r(x)}\left( \int W_r(\xi) P(x|\xi) d\xi \right) dx\\
&=\int \frac{H(x)}{H_r(x)}H_r(x) dx\\
&=\int H(x) dx
\end{aligned}
\tag{12}
となり、真の画像の次の反復の$W_{r+1}$の積分値は、（$W_r$の積分値によらず、）観測画像$H$と等しくなる。
（$r=0$からはじまる式だとすると、初期値の$W_0$以外は、観測画像と等しい積分値であり、そのため、数値計算において、計算量削減のため観測画像の合計によらない$W_0$が入力される場合がある。）
RL法のPSFの規格化性
RL法は、式の特性上、PSFの合計が$1$に規格化される。この性質について、式を追って簡単に説明する。
PSFの合計が$1$の$P$と、その$\alpha$倍したPSFを$P_\alpha$とする。$P_\alpha =\alpha P$となるので、式(11)は、
\begin{aligned}
W_{r+1}(\xi)&=W_r(\xi)\int \frac{H(x)}{H_r(x)}P_\alpha(x|\xi)dx\\
&=W_r(\xi)\int H(x) \frac{P_\alpha(x|\xi)}{\int W_r(\xi)P_\alpha(x|\xi)d\xi}dx\\
&=W_r(\xi)\int H(x) \frac{\alpha P(x|\xi)}{\int W_r(\xi)\alpha P(x|\xi)d\xi}dx\\
&=W_r(\xi)\int H(x) \frac{ P(x|\xi)}{\int W_r(\xi) P(x|\xi)d\xi}dx\\
\end{aligned}
\tag{13}
となる。式(13)から、PSFの合計が$\alpha$の$P_\alpha$は、PSFの合計が$1$の$P$と同じ計算結果となる。しかし、これだけでは、$P$がたまたま$1$であっただけ、すなわち、PSFの合計が任意の規格化が行われているだけで、必ずしも$1$で規格化されているとは限らない。そこで、式(12)と同様に、反復回数毎の推定される真の画像$W_{r+1}$の合計が$H$と等しくなることを確認することで、PSFが$1$に規格化されていることを示す。式(11)の両辺を$P_\alpha$を代入して積分すると、
\begin{aligned}
\int W_{r+1}(\xi)d\xi&=\int \left(W_r(\xi)\int \frac{H(x)}{H_r(x)}P_\alpha(x|\xi)dx \right) d\xi\\
&=\int \frac{H(x)}{H_r(x)}\left( \int W_r(\xi) P_\alpha(x|\xi) d\xi \right) dx\\
&=\int \frac{H(x)}{\int W_r(\xi) P_\alpha(x|\xi)d\xi}\left( \int W_r(\xi) P_\alpha(x|\xi) d\xi \right) dx\\
&=\int \frac{H(x)}{\int W_r(\xi) \alpha P(x|\xi)d\xi}\left( \int W_r(\xi) \alpha P(x|\xi) d\xi \right) dx\\
&=\int H(x) dx
\end{aligned}
\tag{14}
となる。このように、$P_\alpha$は$\alpha$倍したときと、$W_{r+1}$と$H$の合計が一致することから、PSFは1に規格化されている。（例えば、観測に欠落があることをPSFで表現しようとしてPSFの合計が0.5としてRL法を行っても、推定される真の画像は観測画像と等しくなるので、そのような補正はRL法を行った後に別途規格化する必要がある。また、数値計算においては、PSFを合計を1に規格化する必要がないので、この箇所の計算を省略できる。）
RL法の収束性
ここで、RL法の収束性について説明する。$H(x)/H_r(x)\sim 1$のとき、式(11)は、
\begin{aligned}
W_{r+1}(\xi)&=W_r(\xi)\int \frac{H(x)}{H_r(x)}P(x|\xi)dx\\
&\sim W_r(\xi)\int P(x|\xi)dx\\
&\sim W_r(\xi)
\end{aligned}
\tag{15}
となる。式(15)から、$H_r(x)\sim H(x)$のとき、$W_{r+1}\sim W_r$となり、RL法の真の画像は収束する。
RL法を図で整理すると以下のようなイメージである。
このように、RL法は、擬似的な観測画像と本物の観測画像に近づくような真の画像をベイズ推定により求める手法である。そのため、畳み込みで計算で表せていない本物の観測画像の特徴（ノイズ）が入っている場合には、正しく真の画像が推定されないので注意したい。
反復処理の定性的な理解
RL法の反復処理のイメージについて数式と図を使って説明する。式(11)を改めて書き下し、注目するところを太字で書き表すと、
W_{r+1}(\xi)=W_r(\xi)\int \frac{H(x)}{H_r(x)}\mathbf{P(x|\xi)}dx\ \quad r = 0,\ 1,\ 2,\ \dots \tag{16}
となる。$H(x)/H_r(x)$の項は、本物の観測画像と擬似的な観測画像の比較を行っている部分である。ただし、その比をそのまま信頼してはならず、というのも、今比較しているのは観測画像同士であり、本来知りたいのは真の画像の確からしさであるからである。真の画像は、PSFの確率分布で広がって観測されているはずなので、その周辺での$H/H_r$の比も考慮する必要があり、それはPSFの重み($P(x|\xi)$)をかけることで行われ、図で示すと以下のようなイメージで更新を行なうことになる。
※以下のイメージ図は、本質的には同じようなことなので、簡略化のため$H/H_r$ではなく、$H-H_r$をイメージしたものである。
このように、反復初期では、推定される真の画像$W_r$が大きく異なるため、$H/H_r$も非常に異なり、大きく更新されることがわかる。では、反復がある程度進むとどうなるだろうか。以下のようになる。
反復が進むと、実際の観測画像と擬似的に作った観測画像が近づいてくるため、微調整のフェーズに入る。
この段階では、軽微なノイズ等のプラス/マイナスで現れる値などは、$H/H_r$比がPSFで畳み込まれるので、平滑化されて打ち消し合い、真の画像の更新に影響をあまり与えないで収束に向かう。
補足ですが、$H/H_r$比がPSFで畳み込まれるため、より真の画像のパラメータを更新して近づくのは、やはりPSFが小さい場合である。PSFで畳み込む計算は、PSFよりも細かい構造を予想するのは、手法的に困難となることも意味している。
このことについて、簡単にシミュレーションで確かめる。以下がその結果である。
図は、真の画像をPSFの広がりが小さい場合と大きい場合で、畳み込みにより観測画像を用意し、それぞれの画像に対して畳み込んだPSFを使ってRL法を行った結果である。図から、PSFが小さい上の画像の方が推定される真の画像が、実際の真の画像に近いことが確認できる。
このように、RL法はPSFが小さいほど、真の画像に近いものを生成できる可能性が高まることがわかる。究極的には、PSFが単写、すなわちブレが完全にない場合であり、その場合の真の画像の振る舞いも、以下でまとめた。
真の画像初期値依存性
反復処理の定性的な理解で、PSFよりも細かい構造を予想するのは、手法的に困難となることも意味していると述べたが、この事から、真の画像の初期値$W_0$には、PSFより細かい構造を持った分布を入れない方が良いと考えられる。
この真の画像の初期値の依存性について、シミュレーションして確かめる。ここでは、真の画像が点源であり、PSFの広がりがある程度大きい分布に対して、PSFで畳み込んで観測画像を用意し、その画像に対してRL法を行う。下図がその結果である。
- 左図：観測のシミュレーションをするために、PSFで真の画像を畳み込んだイメージ図である。
- 右図：そのPSFを使ってRL法を行った結果で、真の画像の初期値$W_0$にそれぞれ異なる値を入れた画像を使って比較した図である。
このように、$W_0$の初期値にPSFより細かい構造を入れた場合、観測画像での$H/H_r$比に対して、PSFで畳み込むため、PSFの広がり以下の更新の適切な重みが得られず、左図の実際の真の画像と比べて、$W_0$に依存した際立った真の画像が推定されることがわかる。このことから、$W_0$には全てのピクセル値が同じ値や観測画像などの使用が望ましいと考えられる。
PSFの特殊ケース(単写の場合の挙動)
極端なケースとしてPSFが単写、すなわちブレが全くない場合の特殊なケースのRL法の振る舞いについて一応考える。この関係を、PSFはディラックのデルタ関数を用いて表すと、
\begin{aligned}
H(x)&=\int W(\xi)P(x|\xi) d\xi \\
&=\int W(\xi)P(x-\xi) d\xi \\
&= \int W(\xi)\delta (x-\xi ) d\xi \\
&=W(x)
\end{aligned}
\tag{17}
となる。（積分の記号を使っているため、無限に小さい区間$[x、x+dx]$で確率密度が1を満たすようにディラックのデルタ関数を使った。）反復回数$r=0$での真の画像と観測画像で必ず$W_0=H_0$が成り立ち、これをRL法の式(11)に代入して式(17)を考慮して計算すると、
\begin{aligned}
W_{1}(\xi)&=W_0(\xi)\int \frac{H(x)}{H_0(x)}P(x|\xi)dx\\
&=W_0(\xi)\int \frac{H(x)}{W_0(x)}\delta (x-\xi )dx\\
&=W_0(\xi)\frac{H(\xi)}{W_0(\xi)}\\
&=H(\xi)
\end{aligned}
\tag{18}
となる。つまり、PSFが単写の場合、初期状態によらず、1反復で確実に観測画像に収束することになる。このような単写のPSFでもRL法は正しく計算される。
最尤法とRL法の関係
RL法と最尤法の関係について式を追いながら説明する。最尤法とは、データの分布やモデルの仮定に基づいて、データをもっともよく説明するパラメータを推定する方法である。つまり、尤度関数を最大化するようなパラメータを求める方法である。RL法では、観測画像$H$を説明する最も尤もらしい真の画像$W$を求めることに対応する。
尤度関数は、真の画像$W$の各画素の$n$個のピクセル値を$W=[W(0),\ W(1),\ \ldots,\ W(n)]$とし、同じサイズの観測画像$H$の各画素のピクセル値を$H=[H(0),\ H(1),\ \ldots,\ H(n)]$とすると、尤度関数は、
L(W| H) = \prod_{x=1}^{n} P(H(x)|W) \tag{19}
となる。（$\prod$は、直積を表し、例えば$\prod_{a=1}^{3}a_i=a_1\cdot a_2\cdot a_3$という具合である。）
最尤法では、尤度関数$L$を最大化するようなパラメータ$\hat{W}$（最尤推定量と呼ばれる）を求めることが目的であり、
\begin{aligned}
\hat{W} &= \arg\max_{W} L(W| H)\\
&=\arg\max_{W} \prod_{x=1}^{n} P(H(x)|W)
\end{aligned}
\tag{20}
と表せる。（$\arg\max_{W}$とは、ある関数（ここでは、尤度関数$L$）を最大化させるようなパラメータ$W$を求めることを意味する。)
式(20)は直積の計算であり計算が複雑なので、対数(log)を取り計算を簡易化すると、
\begin{aligned}
\hat{W} &= \arg\max_{W} \prod_{x=1}^{n} P(H(x)|W)\\
&= \arg\max_{W} \log \prod_{x=1}^{n} P(H(x)|W)\\
&= \arg\max_{W} \sum_{x=1}^{n} \log P(H(x)|W)
\end{aligned}
\tag{21}
となる。（対数は単調増加関数なので、logを取る前と後で、最小値、最大値などの対応関係は変わらない。）
この式(21)を$W$の各変数毎で偏微分してその傾きが$0$として解いた極値が、適切な仮定のもとでしばしば極大値でかつ最大値となる最尤推定量$\hat{W}$となり、
\nabla \sum_{x=1}^{n} \log P(H(x)|W)=\mathbf{0}\\
\tag{22}
と表される。ここで、$\nabla =\left(\frac {\partial }{\partial W(0)},\ \frac {\partial }{\partial W(1)},\ \dots,\ \frac {\partial }{\partial W(n)}\right)$、n次元の零ベクトルを$\mathbf{0}$とする。この式のある$k$番目に着目して整理すると、
\begin{aligned}
\frac {\partial }{\partial W(k) } \sum_{x=1}^{n} \log P(H(x)|W)&=0\\
\Leftrightarrow \sum_{x=1}^{n} \frac {\partial }{\partial W(k) } \log P(H(x)|W)&=0\\
\Leftrightarrow \int \frac {\partial }{\partial W(k) }\log P(H(x)|W)dx&=0\\
\ \quad (k = 0,\ 1,\ 2,\ \dots,\ n)
\end{aligned}
\tag{23}
と表される。ここで、積分記号は数式を簡易化するために用いた。この方程式は最尤方程式と呼ばれる。
尤度関数がポアソン分布の場合
カウント画像など統計的に観測を得る場合、ポアソンノイズが発生するので、尤度関数はポアソン分布に従い、
L(W| H) = \prod_{x=1}^{n} P(H(x)|W) =\prod_{x=1}^{n} \frac{(W \otimes P)(x)^{H(x)}e^{-(W \otimes P)(x)}}{H(x)!}
\tag{24}
と表される。ここで、$(W \otimes P)(x)$は、真の画像$W$とPSF$P$の畳み込んだ画像（擬似的な観測画像）の$x$番目のピクセル値を表す。
最尤推定量$\hat{W}$を求めるため、真の画像の任意の$k$番目のピクセルについて式(23)の左辺に、式(24)を代入して計算すると、
\begin{aligned}
(左辺)&=\int \frac {\partial }{\partial W(k) } \log P(H(x)|W)dx\\
&=\int \frac {\partial }{\partial W(k) }\log \frac{(W \otimes P)(x)^{H(x)}e^{-(W \otimes P)(x)}}{H(x)!}dx\\
&=\int \frac {\partial }{\partial W(k) }\left[H(x)\log (W \otimes P)(x)-(W \otimes P)(x)-\log H(x)!\right]dx\\
&=\int \frac {\partial }{\partial (W \otimes P)(x)}\left[H(x)\log (W \otimes P)(x)-(W \otimes P)(x)-\log H(x)!\right]\frac {\partial (W \otimes P)(x)}{\partial W(k)}dx\\
&=\int \left( \frac {H(x)}{(W \otimes P)(x)}-1\right)\frac {\partial (W \otimes P)(x)}{\partial W(k)}dx\\
&=\int \left( \frac {H(x)}{(W \otimes P)(x)}-1\right)\frac {\partial}{\partial W(k)}\int W(\xi)P(x-\xi)d\xi dx\ (\because \ 式(4))\\
&=\int \left( \frac {H(x)}{(W \otimes P)(x)}-1\right)P(x-k)dx\\
&=\int \frac {H(x)}{(W \otimes P)(x)}P(x-k)dx-\int P(x-k)dx\\
&=\int \frac {H(x)}{(W \otimes P)(x)}P(x-k)dx-1\ \left(\because \ \int P(x-k)dx=1\right)\\
&=\int F_1(x)P(x-k)dx-1\ \left(F_1(x)\colon =\frac {H(x)}{(W \otimes P)(x)}\right)\\
&=\int F_1(x)P^*(k-x)dx-1\ (P^*はPの各軸に対して反転した配列)\\
&=(F_1 \otimes P^*)(k)-1\\
\end{aligned}
\tag{25}
となる。式(25)は式(23)より
(F_1 \otimes P^*)(k)-1=0
\tag{26}
となる。ここで、反復回での収束条件$W_{r+1}(k)/W_r(k)\sim 1\ (r = 0,\ 1,\ 2,\ \dots)$を導入すると式(26)は、
\begin{aligned}
(F_1\otimes P^*)(k)-\frac{W_{r+1}(k)}{W_r(k)}&=0\\
\Leftrightarrow W_{r+1}(k)&=W_r(k)(F_1\otimes P^*)(k)\\
\Leftrightarrow W_{r+1}(k)&=W_r(k)\left(\frac {H}{W_r \otimes P}\otimes P^*\right)(k)\\
\end{aligned}
\tag{27}
となり、RL法の式(5)が導かれる。このように、RL法は最尤法を実装しており、ポアソンノイズの存在下で良質なデコンボリューションを行う手法である。
尤度関数がガウス分布の場合
尤度関数がガウス分布についても考える。デジタルカメラなどの撮影の場合、観測画像にガウシアンノイズが発生し、尤度関数はガウス分布に従い、
L(W| H) = \prod_{x=1}^{n} P(H(x)|W) =\prod_{x=1}^{n} \frac{1}{\sqrt{2\pi \sigma^2}}\exp \left(-\frac{(H(x)-(W \otimes P)(x))^2}{2\sigma^2}\right)
\tag{28}
と表される。ここで、$(W \otimes P)(x)$は、真の画像$W$とPSF$P$の畳み込んだ画像（擬似的な観測画像）の$x$番目のピクセル値を表す。$\sigma^2$は、分散とする。
最尤推定量$\hat{W}$を求めるため、任意の$k$について式(23)の左辺に、式(28)を代入して計算すると、
\begin{aligned}
(左辺)&=\int \frac {\partial }{\partial W(k) } \log P(H(x)|W)dx\\
&=\int \frac {\partial }{\partial W(k) }\log \left(\frac{1}{\sqrt{2\pi \sigma^2}}\exp \left(-\frac{(H(x)-(W \otimes P)(x))^2}{2\sigma^2}\right)\right) dx\\
&=\int \frac {\partial }{\partial W(k) }\left(\frac{1}{\sqrt{2\pi \sigma^2}}-\frac{(H(x)-(W \otimes P)(x))^2}{2\sigma^2}\right) dx\\
&=\int \frac {\partial }{\partial (W \otimes P)(x)}\left(\frac{1}{\sqrt{2\pi \sigma^2}}-\frac{(H(x)-(W \otimes P)(x))^2}{2\sigma^2}\right) \frac {\partial (W \otimes P)(x)}{\partial W(k)} dx\\
&=\frac{1}{\sigma^2} \int (H(x)-(W \otimes P)(x)) \frac {\partial (W \otimes P)(x)}{\partial W(k)} dx\\
&=\frac{1}{\sigma^2} \int (H(x)-(W \otimes P)(x)) \frac {\partial}{\partial W(k)}\int W(\xi)P(x-\xi)d\xi dx\ (\because \ 式(4))\\
&=\frac{1}{\sigma^2} \int (H(x)-(W \otimes P)(x)) P(x-k) dx\\
&=\frac{1}{\sigma^2} \int (H(x)-(W \otimes P)(x)) P^*(k-x) dx\\
&=\frac{1}{\sigma^2} \int F_2(x) P^*(k-x) dx\ \left(F_2(x)\colon =H(x)-(W \otimes P)(x)\right)\\
&=\frac{1}{\sigma^2} (F_2 \otimes P^* )(k)
\end{aligned}
\tag{29}
となる。式(29)は式(23)より
\begin{aligned}
\frac{1}{\sigma^2}(F_2 \otimes P^* )(k)&=0\\
\Leftrightarrow (F_2 \otimes P^* )(k)&=0
\end{aligned}
\tag{30}
となる。ここで、反復回での収束条件$W_{r+1}(k)-W_r(k)\sim 0\ (r = 0,\ 1,\ 2,\ \dots)$を導入すると式(30)は、
\begin{aligned}
(F_2 \otimes P^* )(k)&=W_{r+1}(k)-W_r(k)\\
\Leftrightarrow W_{r+1}(k)&=W_r(k)+(F_2 \otimes P^* )(k)\\
\Leftrightarrow W_{r+1}(k)&=W_r(k)+\left[ (H-W_r \otimes P) \otimes P^* \right](k)\\
\end{aligned}
\tag{31}
となる。このように、仮定する尤度関数によって得られる式が異なることがわかる。式(31)は、反復回数とポアソン分布の時のRL法の導出過程と同じように解いたものなので、ガウシアンノイズでのRL法として知られている。
おわりに
本記事では、Richardson-Lucy deconvolution（RL法）という古典的な手法について、できるだけ丁寧に解説しました。
最後に、理論的にRL法である程度真の画像が推定できるというイメージが掴めてきたと思うが、実際はもっと複雑で、PSFの場所依存やノイズの影響など様々な問題がある。ここではポアソンノイズについて定性的にその難しさを説明する。
このポアソンノイズを理解するために、私は少しダーツをやったこともあり、ダーツボードを観測画像だと思って説明する。とりあえず、ビギナーの方にどこか一点狙いをしてもらうという設定で、第三者に予想してもらうと以下のようになるだろう。
答えは19のトリプルであるが、このように、ダーツを投げてもらう回数に応じて、どこを狙っているがわかってくる。一方、同じことをプロの方に投げてもらうと以下のようになるだろう。
答えは20のトリプルであるが、数回投げただけでどこを狙っているかわかる。そこで、事前に相手の投げる癖をよく知り正確なPSFがわかっていれば、RL法を使って試合を有利に進められる。というのは冗談だが、他のイメージデコンボリューションでも同じく、PSFが大きいとそれだけ真の画像の候補が増えるため、PSFの大きさと真の画像の推定精度において重要な役割を持つ。さらに、統計誤差をはじめとしたノイズが大きい場合は、観測画像を完全に信頼するのではなく、ノイズとの影響も考慮する必要がある。また、実際のブレ除去問題では、PSFを正確に取得することが難しい場合が多くある。このように、現実のブレ除去は非常に複雑な問題であり、RL法は、あくまでイメージデコンボリューションの古典的な手法として位置付けられ、現在ではさまざまな新しい手法が日々研究されている。
参考文献
[1] William Hadley Richardson. Bayesian-based iterative method of image restoration. JoSA, Vol. 62, No. 1, pp. 55–59, 1972.
[2] Leon B Lucy. An iterative technique for the rectification of observed distributions. The astronomical journal, Vol. 79, p. 745, 1974.
[3] http://hvrl.ics.keio.ac.jp/charmie/doc/ImageRestoration/RichardsonLucy.pdf
[4] http://hvrl.ics.keio.ac.jp/charmie/doc/ImageRestoration/RichardsonLucy_Gaussian.pdf