撮像対象の画像空間を考慮した再構成法

コンプトンカメラは、各測定事象に対応するシステム応答関数を円錐面として取得す る。そのため、コンプトンカメラではPETやSPECTのようにシステム応答関数が直線状 である技術よりも、画像空間とシステム応答関数との交差領域が測定事象によって大きく 異なる。 この現象は(9)および(32)の右辺の分母の大きさに影響を与え、画像空間とシス テム応答関数との交差領域が小さくなるほど画像への寄与を著しく増大させる。

そこで本研究では、LM-ML-EMアルゴリズムに位置依存的に変化する画像空間とコン プトンコーンとの交差領域の影響を組み込むことにより、測定事象によって変化する再構 成画像への寄与を抑制する画像再構成法を開発した。

以下のように、空間中の線源分布𝝀の最尤推定量を得るために、最大化する尤度関数を 定義する:

𝐿(𝔸,𝒥,𝑁|𝑇,𝝀) = ln𝑝(𝑁|𝑇,𝝀) + ln𝑝(𝔸,𝒥|𝝀), (56) 𝐿(𝔸,𝒥,𝑇|𝑁,𝝀) = ln𝑝(𝑇|𝑁,𝝀) + ln𝑝(𝔸,𝒥|𝝀). (57) これらは(17)と(18)に対応する。(17)および(18)との違いは、(56)および(57)では𝔸を得る確 率 [ここでは𝑝(𝔸,𝒥|𝝀)] に撮像対象の画像空間を表す𝒥を含んでいることである。

𝑝(𝔸,𝒥|𝝀)を潜在変数𝕫(20)とKL divergenceを用いて以下のように表す:

40

𝑝(𝔸,𝒥|𝝀) =∏ ∑ 𝑝(𝑧_𝑖𝑗,𝑨_𝑖,𝒥|𝝀)

𝑀 𝑗=1 𝑁 𝑖=1

, (58)

ln𝑝(𝔸,𝒥|𝝀) =ℒ[𝑞(𝕫),𝝀] + KL[𝑞(𝕫)∥ 𝑝(𝕫,𝒥|𝔸,𝝀)], (59) ℒ[𝑞(𝕫),𝝀] =∑ ∑ 𝑞(𝑧_𝑖𝑗)

𝑀 𝑗=1 𝑁 𝑖=1

ln[𝑝(𝑧_𝑖𝑗,𝑨_𝑖,𝒥|𝝀)

𝑞(𝑧_𝑖𝑗) ], (60)

KL[𝑞(𝕫)∥ 𝑝(𝕫,𝒥|𝔸,𝝀)] =− ∑ ∑ 𝑞(𝑧_𝑖𝑗)

𝑀 𝑗=1 𝑁 𝑖=1

ln[𝑝(𝑧_𝑖𝑗,𝒥|𝑨_𝑖,𝝀)

𝑞(𝑧_𝑖𝑗) ]. (61)

∵ 𝑞(𝕫) = {𝑞(𝑧₁₁), … ,𝑞(𝑧_𝑁𝑀)},

𝑝(𝕫,𝒥|𝔸,𝝀) = {𝑝(𝑧₁₁,𝒥|𝑨₁,𝝀), … ,𝑝(𝑧_𝑁𝑀,𝒥|𝑨_𝑁,𝝀)}.

線源分布の𝑙回目の更新値𝝀^(𝑙)を得ている場合、KL = 0となる条件は𝑞(𝑧_𝑖𝑗)=

𝑝(𝑧_𝑖𝑗,𝒥|𝑨_𝑖,𝝀^(𝑙))である。したがって、(57)を例としたとき、E-stepの結果は次式となる: 𝑄(𝝀|𝝀^(𝑙))+ const.

= ln𝑝(𝑇|𝑁,𝝀) +∑ ∑ 𝑝(𝑧_𝑖𝑗,𝒥|𝑨_𝑖,𝝀^(𝑙))ln𝑝(𝑧_𝑖𝑗,𝑨_𝑖,𝒥|𝝀)

𝑀 𝑗=1 𝑁 𝑖=1

− ∑ ∑ 𝑝(𝑧_𝑖𝑗,𝒥|𝑨_𝑖,𝝀^(𝑙))ln𝑝(𝑧_𝑖𝑗,𝒥|𝑨_𝑖,𝝀^(𝑙))

𝑀 𝑗=1 𝑁 𝑖=1

.

(62)

左辺で定数項として表現しているのは、𝝀に依存しない右辺の最後の項である。

M-stepの十分条件は𝑄の𝝀による偏微分導関数が0となる場合であり、以下のように表

す:

∇_𝝀𝑄(𝝀|𝝀^(𝑙))=∇_𝝀ln𝑝(𝑇|𝑁,𝝀) +∇_𝝀∑ ∑ 𝑝(𝑧_𝑖𝑗,𝒥|𝑨_𝑖,𝝀^(𝑙))ln𝑝(𝑧_𝑖𝑗,𝑨_𝑖,𝒥|𝝀)

𝑀 𝑗=1 𝑁 𝑖=1

. (63) 以下、具体的なアルゴリズム導出のために、これまでに現れた確率および変数それぞれ を説明する。

𝑝(𝑧_𝑖𝑗,𝑨_𝑖,𝒥|𝝀)を以下のように表す:

𝑝(𝑧_𝑖𝑗,𝑨_𝑖,𝒥|𝝀)=𝑝(𝑗,𝒥,𝑨_𝑖|𝝀)^𝑧𝑖𝑗. (64) 𝑝(𝑗,𝒥,𝑨_𝑖|𝝀) =𝑝(𝑗|𝒥,𝝀)𝑝(𝑨_𝑖|𝑗)𝑝(𝒥). (65)

𝑝(𝑗,𝒥,𝑨_𝑖|𝝀)は、線源分布𝝀が含まれる画像空間𝒥を測定するときに、得られた光子測定デ

ータ𝑨_𝑖が画素𝑗から放出された光子によって引き起こされた確率である。𝑝(𝑗,𝒥,𝑨_𝑖|𝝀)の成 分である𝑝(𝑗|𝒥,𝝀)は、ある測定事象が𝒥に含まれる画素𝑗から放出される光子によって引き 起こされる確率であり、ここでは次式のようにみなす:

𝑝(𝑗|𝒥,𝝀) =𝑝(𝑗|𝝀). (66)

41

これは(31)より、𝑝(𝑗|𝝀)の正規化は、画像空間内の期待される測定事象全ての総和によっ て行われるため、𝒥はすでに条件として組み込まれているためである。𝑝(𝒥)は、𝒥から到 来する光子が有効な測定事象として得られる確率であり、次式のように表す:

𝑝(𝒥) =∑ 𝑠_𝑗

𝑀 𝑗=1

. (67)

しかし、𝑝(𝒥)は𝝀には依存しないのでM-stepにおいて0となる。

𝑝(𝑧_𝑖𝑗,𝒥|𝑨_𝑖,𝝀^(𝑙))は、𝝀^(𝑙)を測定して得られた測定事象𝑨_𝑖が、𝒥に含まれる𝑗から放出され た光子に起因する確率である。これは𝑧_𝑖𝑗の推定値を表すため𝑧_𝑖𝑗を用いずに、以下のよう に表す:

𝑝(𝑧_𝑖𝑗,𝒥|𝑨_𝑖,𝝀^(𝑙))=𝑝(𝑗|𝑨_𝑖,𝝀^(𝑙))𝑝(𝑗,𝒥|𝑨_𝑖). (68) 𝑝(𝑗|𝑨_𝑖,𝝀^(𝑙))は(32)に示す確率と同様である。𝑝(𝑗,𝒥|𝑨_𝑖)は𝑨_𝑖が𝒥に含まれる𝑗から放出され た光子に起因して得られる確率である。

以上より、E-stepの結果を以下のように表す: 𝑄(𝝀|𝝀^(𝑙))= ln𝑝(𝑇|𝑁,𝝀)

+∑ ∑ 𝑝(𝑧_𝑖𝑗,𝒥|𝑨_𝑖,𝝀^(𝑙))[ln𝑝(𝑗|𝒥,𝝀) + ln𝑝(𝑨_𝑖|𝑗) + ln𝑝(𝒥)]

𝑀 𝑗=1 𝑁 𝑖=1

= ln𝑁

𝑇 +𝑁 ln(𝑇 ∑ 𝑠_𝑗𝜆_𝑗

𝑀 𝑗=1

) − 𝑇 ∑ 𝑠_𝑗𝜆_𝑗

𝑀 𝑗=1

−ln𝑁!

+∑ ∑ 𝑝(𝑗|𝑨_𝑖,𝝀^(𝑙))𝑝(𝑗,𝒥|𝑨_𝑖)[ln 𝑠_𝑗𝜆_𝑗

∑^𝑀 𝑠_𝑘𝜆_𝑘

𝑘=1

+ ln𝑝(𝑨_𝑖|𝑗) + ln𝑝(𝒥)]

𝑀 𝑗=1 𝑁 𝑖=1

.

(69)

したがって、M-stepは次式となる:

∇_𝝀𝑄(𝝀|𝝀^(𝑙))=∇_𝝀𝑁 (ln𝑇 + ln∑ 𝑠_𝑗𝜆_𝑗

𝑀 𝑗=1

) − ∇_𝝀𝑇 ∑ 𝑠_𝑗𝜆_𝑗

𝑀 𝑗=1

+∇_𝝀∑ ∑ 𝑝(𝑗|𝑨_𝑖,𝝀^(𝑙))𝑝(𝑗,𝒥|𝑨_𝑖)(ln𝑠_𝑗𝜆_𝑗−ln∑ 𝑠_𝑘𝜆_𝑘

𝑀 𝑘=1

)

𝑀 𝑗=1 𝑁 𝑖=1

=−𝑇 𝑠_𝑗+ 1

𝜆_𝑗∑ 𝑝(𝑗|𝑨_𝑖,𝝀^(𝑙))𝑝(𝑗,𝒥|𝑨_𝑖)

𝑁 𝑖=1

= 0.

(70)

上式の𝜆_𝑗を𝜆_𝑗^(𝑙+1)とすることで、以下の再構成アルゴリズムを得る:

𝜆_𝑗^(𝑙+1) =𝜆_𝑗^(𝑙)

𝑠_𝑗 ∑ 𝑣_𝑖𝑗𝑡_𝑖𝑗

∑^𝑀 𝑡_𝑖𝑘𝜆_𝑘^(𝑙)

𝑘=1 𝑁

𝑖=1

, ∵ 𝑣_𝑖𝑗 =𝑝(𝑗,𝒥|𝑨_𝑖), 𝑡_𝑖𝑗 =𝑝(𝑗,𝑨_𝑖). (71)

42

ここでは(55)と同様に測定の時間依存性を無視している。𝑣_𝑖𝑗の具体的な導出は2.3節にて 説明する。