Procedure for recovering a blurred image
インパルス応答関数によってぼやけた信号を回復するためにリチャードソン・ルーシー逆畳み込みを使用します。 リチャードソン・ルーシーアルゴリズム（ ルーシー・リチャードソン・デコンボリューション とも呼ばれる）は、 既知の 点像分布関数によって ぼかさ れた画像を復元するための 反復的な手順 である。 このアルゴリズムは、独立してこのアルゴリズムを開発したウィリアム・リチャードソンと レオン・B・ルーシー にちなんで名付けられた 。 [1] [2]
説明 光学システムを用いて画像が生成され、 写真フィルム 、 電荷結合素子、 CMOSセンサー などを用いて検出される場合 、画像は必然的にぼやけ、理想的な 点光源 は点として現れず、点像分布関数と呼ばれる広がり方をします。広がりのある光源は多数の点光源の総和に分解できるため、観測画像は、基となる画像に作用する
遷移行列 p として表すことができます。
d
i
=
∑
j
p
i
,
j
u
j
{\displaystyle d_{i}=\sum _{j}p_{i,j}u_{j}\,}
ここで 、 はピクセル における基礎画像の強度 、 は ピクセル で検出された強度です 。一般に、 を要素とする行列は、ピクセル i で検出された光源ピクセル j からの光の量を表します。ほとんどの優れた光学システム（または一般に、 シフト不変 と記述される線形システム ）では、伝達関数 p は 光源ピクセル j と観測ピクセル i 間の
空間 オフセット で単純に表すことができます。
u
j
{\displaystyle u_{j}}
j
{\displaystyle j}
d
i
{\displaystyle d_{i}}
i
{\displaystyle i}
p
i
,
j
{\displaystyle p_{i,j}}
p
i
,
j
=
P
(
i
−
j
)
{\displaystyle p_{i,j}=P(i-j)}
ここでは 点像関数 と呼ばれます 。その場合、上式は 畳み込み になります。これは1次元空間を想定して書かれていますが、ほとんどのイメージングシステムは2次元であり、光源、検出画像、点像関数はいずれも2つのインデックスを持ちます。したがって、2次元検出画像は、基となる画像と2次元点像関数、そして追加された検出ノイズとの畳み込みです 。
P
(
Δ
i
)
{\displaystyle P(\Delta i)}
P
(
Δ
x
,
Δ
y
)
{\displaystyle P(\Delta x,\Delta y)}
観測された と既知の を推定するために、 反復回数 t に対するの 推定値 ( と呼ばれる) が次のように更新される
反復手順が採用されます。
u
j
{\displaystyle u_{j}}
d
i
{\displaystyle d_{i}}
P
(
Δ
i
x
,
Δ
j
y
)
{\displaystyle P(\Delta i_{x},\Delta j_{y})}
u
j
{\displaystyle u_{j}}
u
^
j
(
t
)
{\displaystyle {\hat {u}}_{j}^{(t)}}
u
^
j
(
t
+
1
)
=
u
^
j
(
t
)
∑
i
d
i
c
i
p
i
j
{\displaystyle {\hat {u}}_{j}^{(t+1)}={\hat {u}}_{j}^{(t)}\sum _{i}{\frac {d_{i}}{c_{i}}}p_{ij}}
どこ
c
i
=
∑
j
p
i
j
u
^
j
(
t
)
{\displaystyle c_{i}=\sum _{j}p_{ij}{\hat {u}}_{j}^{(t)}}
と 仮定する。この反復が収束する場合、最大尤度解に収束することが経験的に示されている 。 [3]
∑
j
p
i
j
=
1
{\displaystyle \sum _{j}p_{ij}=1}
u
j
{\displaystyle u_{j}}
これを、点広がり関数 P との
畳み込み の観点から、2 次元 (またはそれ以上) についてより一般的に記述すると次のようになります。
u
^
(
t
+
1
)
=
u
^
(
t
)
⋅
(
d
u
^
(
t
)
⊗
P
⊗
P
∗
)
{\displaystyle {\hat {u}}^{(t+1)}={\hat {u}}^{(t)}\cdot \left({\frac {d}{{\hat {u}}^{(t)}\otimes P}}\otimes P^{*}\right)}
ここで、除算と乗算は要素ごとに行われ、 2D 畳み込みを示し、 は鏡像 点像分布関数、または光伝達 関数の エルミート転置 の 逆フーリエ変換です 。
⊗
{\displaystyle \otimes }
P
∗
{\displaystyle P^{*}}
点広がり関数が 事前に 分からない問題では、 ブラインドデコンボリューション を実現するためにリチャードソン・ルーシーアルゴリズムの修正版が提案されている 。 [4]
p
i
j
{\displaystyle p_{ij}}
導出 蛍光顕微鏡の文脈では、ピクセルを持つ検出器の 期待値に対する 光子数（または検出された光に比例するデジタル化カウント）のセットを測定する確率は 次のように与えられる。
m
=
[
m
0
,
.
.
.
,
m
K
]
{\displaystyle \mathbf {m} =[m_{0},...,m_{K}]}
E
=
[
E
0
,
.
.
.
,
E
K
]
{\displaystyle \mathbf {E} =[E_{0},...,E_{K}]}
K
+
1
{\displaystyle K+1}
P
(
m
|
E
)
=
∏
i
K
P
o
i
s
s
o
n
(
E
i
)
=
∏
i
K
E
i
m
i
e
−
E
i
m
i
!
{\displaystyle P(\mathbf {m} \vert \mathbf {E} )=\prod _{i}^{K}\mathrm {Poisson} (E_{i})=\prod _{i}^{K}{\frac {{E_{i}}^{m_{i}}e^{-E_{i}}}{m_{i}!}}}
最大尤度推定の文脈では、 絶対値を考慮せずに
尤度関数 の 最大値 を見つけることが目的なので、これは扱いやすいです。
ln
(
P
)
{\displaystyle \ln(P)}
ln
(
P
(
m
|
E
)
)
=
∑
i
K
[
(
m
i
ln
E
i
−
E
i
)
−
ln
(
m
i
!
)
]
{\displaystyle \ln(P(\mathbf {m} \vert \mathbf {E} ))=\sum _{i}^{K}\left[(m_{i}\ln E_{i}-E_{i})-\ln(m_{i}!)\right]}
は定数なので 、最大値の位置に関する追加情報は得られません。そこで、
ln
(
m
i
!
)
{\displaystyle \ln(m_{i}!)}
α
(
m
|
E
)
=
∑
i
K
[
m
i
ln
E
i
−
E
i
]
{\displaystyle \alpha (\mathbf {m} \vert \mathbf {E} )=\sum _{i}^{K}\left[m_{i}\ln E_{i}-E_{i}\right]}
ここで 、 は と同じ最大位置を共有するものである 。ここで、 が グラウンドトゥルース と線形であると仮定された
測定 から得られるとしよう。すると、
α
{\displaystyle \alpha }
P
(
m
|
E
)
{\displaystyle P(\mathbf {m} \vert \mathbf {E} )}
E
{\displaystyle \mathbf {E} }
x
{\displaystyle \mathbf {x} }
H
{\displaystyle \mathbf {H} }
E
=
H
x
{\displaystyle \mathbf {E} =\mathbf {H} \mathbf {x} }
ここで行列の乗算が暗黙的に表される。これは次のようにも書ける。
E
m
=
∑
n
K
H
m
n
x
n
{\displaystyle E_{m}=\sum _{n}^{K}H_{mn}x_{n}}
ここで、真実が どのように混合またはぼかされるかが わかります 。
H
{\displaystyle H}
また、の元の 他の元に対する微分は 次のように書けること
も示せます。
E
{\displaystyle \mathbf {E} }
(
E
i
)
{\displaystyle (E_{i})}
x
j
{\displaystyle x_{j}}
∂
E
i
∂
x
j
=
H
i
j
{\displaystyle {\frac {\partial E_{i}}{\partial x_{j}}}=H_{ij}}
1
ヒント：例えば (5 x 5) の行列と 5 要素の配列 2 つと を 書い て 確認すれば、簡単に理解できます。この最後の式は、 の 1 つの 要素 （例えば 要素）が 他の 要素にどれだけ 影響を与えるかを表す式です （もちろん、 の場合 も考慮されます）。例えば、典型的なケースでは、グラウンドトゥルース 1 つの要素は の近くの要素に影響を与えます が、非常に離れた要素には影響を与えません（これらの行列要素では の値が 期待されます）。
H
{\displaystyle \mathbf {H} }
E
{\displaystyle \mathbf {E} }
x
{\displaystyle \mathbf {x} }
x
{\displaystyle \mathbf {x} }
i
{\displaystyle i}
j
≠
i
{\displaystyle j\neq i}
i
=
j
{\displaystyle i=j}
x
{\displaystyle \mathbf {x} }
E
{\displaystyle \mathbf {E} }
0
{\displaystyle 0}
さて、鍵となる任意のステップは、 知られていないが、によって推定される可能性がある。RL アルゴリズムを使用する際に推定されたグラウンドトゥルース を と 呼びましょう。 ハット シンボルはグラウンドトゥルースとグラウンドトゥルースの推定値を区別するために使用されます。
x
{\displaystyle \mathbf {x} }
x
^
{\displaystyle {\hat {\mathbf {x} }}}
x
^
o
l
d
{\displaystyle {\hat {\mathbf {x} }}_{old}}
x
^
n
e
w
{\displaystyle {\hat {\mathbf {x} }}_{new}}
x
^
n
e
w
=
x
^
o
l
d
+
λ
∂
α
(
m
|
E
(
x
)
)
∂
x
|
x
^
o
l
d
{\displaystyle {\hat {\mathbf {x} }}_{new}={\hat {\mathbf {x} }}_{old}+\lambda {\frac {\partial \ \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}{\partial \mathbf {x} }}|_{{\hat {\mathbf {x} }}_{old}}}
2
ここで、は次元勾配 を表す 。偏微分を行うと、 次の式が得られる
。
∂
∂
x
{\displaystyle {\frac {\partial }{\partial \mathbf {x} }}}
K
{\displaystyle K}
α
(
m
|
E
(
x
)
)
{\displaystyle \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}
∂
α
(
m
|
E
(
x
)
)
∂
x
j
=
∂
∂
x
j
∑
i
K
[
m
i
ln
E
i
−
E
i
]
=
∑
i
K
[
m
i
E
i
∂
∂
x
j
E
i
−
∂
∂
x
j
E
i
]
=
∑
i
K
∂
E
i
∂
x
j
[
m
i
E
i
−
1
]
{\displaystyle {\frac {\partial \ \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\sum _{i}^{K}\left[m_{i}\ln E_{i}-E_{i}\right]=\sum _{i}^{K}\left[{\frac {m_{i}}{E_{i}}}{\frac {\partial }{\partial x_{j}}}E_{i}-{\frac {\partial }{\partial x_{j}}}E_{i}\right]=\sum _{i}^{K}{\frac {\partial E_{i}}{\partial x_{j}}}\left[{\frac {m_{i}}{E_{i}}}-1\right]}
（ １ ）を代入すると、
∂
α
(
m
|
E
(
x
)
)
∂
x
j
=
∑
i
K
H
i
j
[
m
i
E
i
−
1
]
{\displaystyle {\frac {\partial \ \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}{\partial x_{j}}}=\sum _{i}^{K}H_{ij}\left[{\frac {m_{i}}{E_{i}}}-1\right]}
行列の定義により転置行列となることに注意してください 。したがって
H
j
i
T
=
H
i
j
{\displaystyle {H}_{ji}^{T}=H_{ij}}
∂
α
(
m
|
E
(
x
)
)
∂
x
j
=
∑
i
K
H
j
i
T
[
m
i
E
i
−
1
]
{\displaystyle {\frac {\partial \ \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}{\partial x_{j}}}=\sum _{i}^{K}{H}_{ji}^{T}\left[{\frac {m_{i}}{E_{i}}}-1\right]}
3
この方程式は から までのすべての要素にわたって成り立つので 、これらの 方程式は単一のベクトル方程式として簡潔に書き直すことができる。
j
{\displaystyle j}
1
{\displaystyle 1}
K
{\displaystyle K}
K
{\displaystyle K}
∂
α
(
m
|
E
(
x
)
)
∂
x
=
H
T
[
m
E
−
1
]
{\displaystyle {\frac {\partial \ \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}{\partial \mathbf {x} }}={\mathbf {H} ^{T}}\left[{\frac {\mathbf {m} }{\mathbf {E} }}-\mathbf {1} \right]}
ここで 、は行列 、、、 は ベクトルです。ここで一見恣意的に見えますが重要なステップとして、
H
T
{\displaystyle \mathbf {H} ^{T}}
m
{\displaystyle \mathbf {m} }
E
{\displaystyle \mathbf {E} }
1
{\displaystyle \mathbf {1} }
λ
=
x
^
o
l
d
H
T
1
{\displaystyle \lambda ={\frac {{\hat {\mathbf {x} }}_{old}}{\mathbf {H} ^{T}\mathbf {1} }}}
4
ここで は の大きさのベクトル（ 、 、 と 同じ ）であり、除算は要素ごとに行われる。( 3 )と( 4 )を用いると、( 2 )は次のように書き直される。
1
{\displaystyle \mathbf {1} }
K
{\displaystyle K}
m
{\displaystyle \mathbf {m} }
E
{\displaystyle \mathbf {E} }
x
{\displaystyle \mathbf {x} }
x
^
n
e
w
=
x
^
o
l
d
+
λ
∂
α
(
m
|
E
(
x
)
)
∂
x
=
x
^
o
l
d
+
x
^
o
l
d
H
T
1
H
T
[
m
E
−
1
]
=
x
^
o
l
d
+
x
^
o
l
d
H
T
1
H
T
m
E
−
x
^
o
l
d
{\displaystyle {\hat {\mathbf {x} }}_{new}={\hat {\mathbf {x} }}_{old}+\lambda {\frac {\partial \alpha (\mathbf {m} \vert \mathbf {E} (x))}{\partial x}}={\hat {\mathbf {x} }}_{old}+{\frac {{\hat {\mathbf {x} }}_{old}}{{\mathbf {H} ^{T}}\mathbf {1} }}{\mathbf {H} ^{T}}\left[{\frac {\mathbf {m} }{\mathbf {E} }}-\mathbf {1} \right]={\hat {\mathbf {x} }}_{old}+{\frac {{\hat {\mathbf {x} }}_{old}}{\mathbf {H} ^{T}\mathbf {1} }}\mathbf {H} ^{T}{\frac {\mathbf {m} }{\mathbf {E} }}-{\hat {\mathbf {x} }}_{old}}
その結果
x
^
n
e
w
=
x
^
o
l
d
H
T
(
m
E
)
/
H
T
1
{\displaystyle {\hat {\mathbf {x} }}_{new}={\hat {\mathbf {x} }}_{old}\mathbf {H} ^{T}\left({\frac {\mathbf {m} }{\mathbf {E} }}\right)/\mathbf {H} ^{T}\mathbf {1} }
5
ここで、除算は行列の要素ごとの除算を指し、 行列として動作しますが、除算と積（ の後に暗黙的に ）は要素ごとに行われます。また、 は、 と仮定しているため、計算できます。
H
T
{\displaystyle \mathbf {H} ^{T}}
x
^
o
l
d
{\displaystyle {\hat {\mathbf {x} }}_{old}}
E
≡
E
(
x
^
o
l
d
)
=
H
x
^
o
l
d
{\displaystyle \mathbf {E} \equiv \mathbf {E} ({\hat {\mathbf {x} }}_{old})=\mathbf {H} {\hat {x}}_{old}}
- 初期推定値 は既知である（通常は実験データに設定される）
x
^
0
{\displaystyle {\hat {\mathbf {x} }}_{0}}
- 測定 関数 は既知である
H
{\displaystyle \mathbf {H} }
一方、 は実験データです。したがって、式( 5 )を連続的に適用すると、尤度ランドスケープ で上昇する（尤度の勾配の方向に移動するため）ことによって グラウンドトゥルースを推定するアルゴリズムが得られます 。 この導出では、収束することが示されておらず、初期の選択への依存性が示されていません [ 引用が必要 ] 。 式( 2 )は、尤度が増加する方向をたどる方法を提供しますが、対数微分の選択は任意であることに注意してください。 一方、式( 4 ) は、反復の前のステップからの動きに 重み付けする方法を導入します。 この項が( 5 )に存在しない場合、 の場合でも、アルゴリズムは推定で動きを出力することに注意してください 。 ここで使用される唯一の戦略は、尤度を最大化することであるため、画像にアーティファクトが導入される可能性があることに注意する価値があります。 この導出では、グラウンドトゥルースの形状に関する事前知識は 使用されないことに注意してください。
m
{\displaystyle \mathbf {m} }
x
n
e
w
{\displaystyle \mathbf {x} _{new}}
m
=
E
(
x
^
o
l
d
)
{\displaystyle \mathbf {m} =\mathbf {E} ({\hat {\mathbf {x} }}_{old})}
x
{\displaystyle \mathbf {x} }
ソフトウェア
参照
参考文献
^ リチャードソン, ウィリアム・ハドリー (1972). 「ベイズに基づく反復法による画像復元」. アメリカ光学会誌 . 62 (1): 55– 59. Bibcode :1972JOSA...62...55R. doi :10.1364/JOSA.62.000055. ^ Lucy, LB (1974). 「観測分布の整流化のための反復法」. 天文学ジャーナル . 79 (6): 745– 754. Bibcode :1974AJ.....79..745L. doi :10.1086/111605. ^ Shepp, LA; Vardi, Y. (1982)、「エミッショントモグラフィーにおける最大尤度再構成」、 IEEE Transactions on Medical Imaging 、 1 (2): 113– 22、 doi :10.1109/TMI.1982.4307558、 PMID 18238264 ^ Fish DA; Brinicombe AM; Pike ER; Walker JG (1995), "Blind deconvolution by means of the Richardson–Lucy algorithm" (PDF) , Journal of the Optical Society of America A , 12 (1): 58– 65, Bibcode :1995JOSAA..12...58F, doi :10.1364/JOSAA.12.000058, S2CID 42733042, 2019年1月10日時点のオリジナル (PDF)からアーカイブ
![{\displaystyle \mathbf {m} =[m_{0},...,m_{K}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f973d7066ac0ff2e9cc11e67b391e45969cbd064)
![{\displaystyle \mathbf {E} =[E_{0},...,E_{K}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5386eb0a83c9edb6ef9a2c1d9fd923f0e9277f01)
![{\displaystyle \ln(P(\mathbf {m} \vert \mathbf {E} ))=\sum _{i}^{K}\left[(m_{i}\ln E_{i}-E_{i})-\ln(m_{i}!)\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4bfa55b24935ebdafdc4f357c12442c9dffca9a2)
![{\displaystyle \alpha (\mathbf {m} \vert \mathbf {E} )=\sum _{i}^{K}\left[m_{i}\ln E_{i}-E_{i}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cfd88c48e28c88526bf8f5f6a3980dc5a6dc6901)
![{\displaystyle {\frac {\partial \ \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\sum _{i}^{K}\left[m_{i}\ln E_{i}-E_{i}\right]=\sum _{i}^{K}\left[{\frac {m_{i}}{E_{i}}}{\frac {\partial }{\partial x_{j}}}E_{i}\right]=\sum _{i}^{K}{\frac {\partial E_{i}}{\partial x_{j}}}\left[{\frac {m_{i}}{E_{i}}}-1\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ee3d48f1b65897a42c7d6da26f2b7a7d91257b95)
![{\displaystyle {\frac {\partial \ \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}{\partial x_{j}}}=\sum _{i}^{K}H_{ij}\left[{\frac {m_{i}}{E_{i}}}-1\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e5d651929da3a2f4c6d2d6bfdd4f8e1fb0fc4842)
![{\displaystyle {\frac {\partial \ \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}{\partial x_{j}}}=\sum _{i}^{K}{H}_{ji}^{T}\left[{\frac {m_{i}}{E_{i}}}-1\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dde50e7b4f0494e2816ffe6409f3c130b7f439e6)
![{\displaystyle {\frac {\partial \ \alpha (\mathbf {m} \vert \mathbf {E} (\mathbf {x} ))}{\partial \mathbf {x} }}={\mathbf {H} ^{T}}\left[{\frac {\mathbf {m} }{\mathbf {E} }}-\mathbf {1} \right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d0f03b030c484c6f62b9e22e916e9b34233bfc60)
![{\displaystyle {\hat {\mathbf {x} }}_{new}={\hat {\mathbf {x} }}_{old}+\lambda {\frac {\partial \alpha (\mathbf {m} \vert \mathbf {E} (x))}{\partial x}}={\hat {\mathbf {x} }}_{old}+{\frac {{\hat {\mathbf {x} }}_{old}}{{\mathbf {H} ^{T}}\mathbf {1} }}{\mathbf {H} ^{T}}\left[{\frac {\mathbf {m} }{\mathbf {E} }}-\mathbf {1} \right]={\hat {\mathbf {x} }}_{old}+{\frac {{\hat {\mathbf {x} }}_{old}}{\mathbf {H} ^{T}\mathbf {1} }}\mathbf {H} ^{T}{\frac {\mathbf {m} }{\mathbf {E} }}-{\hat {\mathbf {x} }}_{old}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/75f026dc1cc3e16d42accba0ec64debc6ed66b2e)
