JP2024052415A - 核医学診断装置、データ処理方法及びプログラム - Google Patents

Abstract

【課題】画質を向上させること。【解決手段】実施形態に係る核医学診断装置は、計測部と、特定部とを備える。計測部は、第１の検出器で検出される第１のイベントに関する第１の光子数情報と、前記第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを取得する。特定部は、前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１のイベントと前記第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する。【選択図】図１

Description

本明細書及び図面に開示の実施形態は、核医学診断装置、データ処理方法及びプログラムに関する。

ＰＥＴ画像の再構成を行う際において、ＴｏＦ（Ｔｉｍｅ Ｏｆ Ｆｌｉｇｈｔ）カーネルを用いて画像再構成が行われる。ここで、画像再構成に用いられるＴｏＦカーネルの時間幅が、検出器のＴｏＦ性能に整合するものでない場合、再構成後の画像の画質が低下する場合がある。従って、検出器のＴｏＦ性能と、画像再構成に用いるＴｏＦカーネルとが、できるだけ一致することが望ましい。

しかしながら、例えば、同時計数イベント数が少ない場合等では、観測揺らぎが大きくなり、ＴｏＦカーネルの算出に誤差が生じる場合がある。

Eric Berg and Simon R Cherry, "Using convolutional neural networks to estimate time-of-flight from PET detector waveforms", vol. 63, no 2, Phys. Med. Biol, 2018年

本明細書及び図面の開示の実施形態が解決しようとする課題の一つは、画質を向上させることである。ただし、本明細書及び図面に開示の実施形態により解決しようとする課題は上記課題に限られない。後述する実施形態に示す各構成による各効果に対応する課題を他の課題として位置づけることもできる。

実施形態に係る核医学診断装置は、計測部と、特定部とを備える。計測部は、第１の検出器で検出される第１のイベントに関する第１の光子数情報と、前記第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを取得する。特定部は、前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１のイベントと前記第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する。

図１は、実施形態に係る核医学診断装置の構成の一例を示した図である。 図２は、実施形態に係る核医学診断装置が行う処理の流れを示したフローチャートである。 図３は、図２のステップＳ３０の処理の一例をより詳細に説明したフローチャートである。 図４は、実施形態に係る核医学診断装置が行う処理について説明した図である。 図５は、図３のステップＳ３００の処理の一例をより詳細に説明したフローチャートである。 図６は、実施形態に係る核医学診断装置が行う処理について説明した図である。 図７は、実施形態に係る核医学診断装置が行う処理について説明した図である。 図８は、実施形態に係る核医学診断装置が行う処理について説明した図である。 図９は、図２のステップＳ３００の処理の別の一例をより詳細に説明したフローチャートである。 図１０は、実施形態に係る核医学診断装置が行う処理について説明した図である。 図１１は、実施形態に係る方法の検証結果の一例について示した図である。 図１２は、実施形態に係る方法の検証結果の一例について示した図である。 図１３は、実施形態に係る方法の検証結果の一例について示した図である。

（実施形態）
以下、図面を参照しながら、核医学診断装置、データ処理方法及びプログラムの実施形態について詳細に説明する。

図１は、実施形態に係る核医学診断装置としてのＰＥＴ装置１００の構成を示す図である。図１に示すように、実施形態に係るＰＥＴ装置１００は、架台装置１と、コンソール装置２とを備える。架台装置１は、検出器３と、フロントエンド回路１０２と、天板１０３と、寝台１０４と、寝台駆動部１０６とを備える。

検出器３は、被検体に投与され集積した薬剤から放出される陽電子が周囲の組織の電子と対消滅を起こすことで発生したガンマ線が発光体と相互作用することにより励起状態となった物質が再び基底状態に遷移する際に再放出される光であるシンチレーション光（蛍光）を検出することにより、放射線を検出する検出器である。また、実施形態においては、検出器３は、チェレンコフ光も検出することができる。検出器３は、被検体に投与され集積した薬剤から放出される陽電子が周囲の組織の電子と対消滅を起こすことで発生したガンマ線の放射線のエネルギー情報を検出する。検出器３は、被写体Ｐの周囲をリング状に取り囲むように複数配置され、例えば複数の検出器ブロックからなる。

検出器３は、典型的には、シンチレータ結晶と、光検出素子からなる光検出面とからなる。

シンチレータ結晶の材質としては、例えばチェレンコフ光を発生させるのに適した材質、例えば、ゲルマニウム酸ビスマス（ＢＧＯ、Ｂｉｓｍｕｔｈ Ｇｅｒｍａｎｉｕｍ Ｏｘｉｄｅ）や、鉛ガラス（ＳｉＯ₂＋ＰｂＯ）、フッ化鉛（ＰｂＦ₂）、ＰＷＯ（ＰｂＷＯ₄）等の鉛化合物を用いることができる。また、別の例として、例えば、ＬＹＳＯ（Ｌｕｔｅｔｉｕｍ Ｙｔｔｒｉｕｍ Ｏｘｙｏｒｔｈｏｓｉｌｉｃａｔｅ）、ＬＳＯ（Ｌｕｔｅｔｉｕｍ Ｏｘｙｏｒｔｈｏｓｉｌｉｃａｔｅ）、ＬＧＳＯ（Ｌｕｔｅｔｉｕｍ Ｇａｄｏｌｉｎｉｕｍ Ｏｘｙｏｒｔｈｏｓｉｌｉｃａｔｅ）等やＢＧＯ等の、シンチレータ結晶を用いてもよい。光検出面３ｂを構成する光検出素子は、例えば複数の画素からなり、これらの画素それぞれは、例えばＳＰＡＤ（Ｓｉｎｇｌｅ Ｐｈｏｔｏｎ Ａｖａｌａｎｃｈｅ Ｄｉｏｄｅ）で構成される。また、検出器３の構成は、上述の例に限られず、一例として、光検出素子は、例えばＳｉＰＭ（Ｓｉｌｉｃｏｎ ｐｈｏｔｏｍｕｌｔｉｐｌｉｅｒ）や、光電子増倍管を用いてもよい。

なお、シンチレータ結晶はモノリシック結晶であってもよく、また光検出素子からなる光検出面は、シンチレータ結晶の６面に例えば配置されてもよい。なお、以下の実施形態では、検出器３のシンチレータ結晶は、モノリシック結晶でない場合をまずは例にとり説明する。

また、架台装置１は、フロントエンド回路１０２により、検出器３の出力信号から計数情報を生成し、生成した計数情報を、コンソール装置２の記憶部１３０に格納する。なお、検出器３は、複数のブロックに区分けされ、フロントエンド回路１０２を備える。

フロントエンド回路１０２は、検出器３の出力信号をデジタルデータに変換し、計数情報を生成する。この計数情報には、消滅ガンマ線の検出位置、エネルギー値、及び検出時間が含まれる。例えば、フロントエンド回路１０２は、シンチレーション光を同じタイミングで電気信号に変換した複数の光検出素子を特定する。そして、フロントエンド回路１０２は、消滅ガンマ線が入射したシンチレータの位置を示すシンチレータ番号（Ｐ）を特定する。消滅ガンマ線が入射したシンチレータ位置を特定する手段は、各光検出素子の位置及び電気信号の強度に基づいて重心演算を行うことによって特定してもよい。また、シンチレータと光検出素子の各々の素子サイズが対応している場合には、例えば最大出力が得られた光検出素子に対応するシンチレータを消滅ガンマ線が入射したシンチレータ位置と仮定し、さらにシンチレータ間散乱を考慮して最終的に特定するなどする。

また、フロントエンド回路１０２は、各光検出素子から出力された電気信号の強度を積分計算するあるいは電気信号強度が閾値を超えた時間（Ｔｉｍｅ ｏｖｅｒ Ｔｈｒｅｓｈｏｌｄ）を計測し検出器３に入射した消滅ガンマ線のエネルギー値（Ｅ）を特定する。また、フロントエンド回路１０２は、検出器３によって消滅ガンマ線によるシンチレーション光が検出された検出時間（Ｔ）を特定する。なお、検出時間（Ｔ）は、絶対時刻であってもよいし、撮影開始時点からの経過時間であってもよい。このように、フロントエンド回路１０２は、シンチレータ番号（Ｐ）、エネルギー値（Ｅ）、及び検出時間（Ｔ）を含む計数情報を生成する。

なお、フロントエンド回路１０２は、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、ＧＰＵ（Ｇｒａｐｈｉｃａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）或いは、特定用途向け集積回路（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ：ＡＳＩＣ）、プログラマブル論理デバイス（例えば、単純プログラマブル論理デバイス（Ｓｉｍｐｌｅ Ｐｒｏｇｒａｍｍａｂｌｅ Ｌｏｇｉｃ Ｄｅｖｉｃｅ：ＳＰＬＤ）、複合プログラマブル論理デバイス（Ｃｏｍｐｌｅｘ Ｐｒｏｇｒａｍｍａｂｌｅ Ｌｏｇｉｃ Ｄｅｖｉｃｅ：ＣＰＬＤ）、及びフィールドプログラマブルゲートアレイ（Ｆｉｅｌｄ Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ：ＦＰＧＡ））等の回路により実現される。フロントエンド回路１０２は、フロントエンド部の一例である。

天板１０３は、被写体Ｐが載置されるベッドであり、寝台１０４の上に配置される。寝台駆動部１０６は、処理回路１５０の制御機能１０５ｅによる制御の下、天板１０３を移動させる。例えば、寝台駆動部１０６は、天板１０３を移動させることで、被写体Ｐを架台装置１の撮影口内に移動させる。

コンソール装置２は、操作者によるＰＥＴ装置１００の操作を受け付け、ＰＥＴ画像の撮影を制御するとともに、架台装置１によって収集された計数情報を用いてＰＥＴ画像を再構成する。図１に示すように、コンソール装置２は、処理回路１５０と、入力装置１１０と、ディスプレイ１２０と、記憶部１３０とを備える。なお、コンソール装置２が備える各部は、バスを介して接続される。処理回路１５０の詳細については後述する。

入力装置１１０は、ＰＥＴ装置１００の操作者によって各種指示や各種設定の入力に用いられるマウスやキーボード等であり、入力された各種指示や各種設定を、処理回路１５０に転送する。例えば、入力装置１１０は、撮影開始指示の入力に用いられる。

ディスプレイ１２０は、操作者によって参照されるモニター等であり、処理回路１５０による制御の下、被写体の呼吸波形やＰＥＴ画像を表示したり、操作者から各種指示や各種設定を受け付けるためのＧＵＩ（Ｇｒａｐｈｉｃａｌ Ｕｓｅｒ Ｉｎｔｅｒｆａｃｅ）を表示したりする。

記憶部１３０は、ＰＥＴ装置１００において用いられる各種データを記憶する。記憶部１３０は、例えば、メモリで構成され、一例として、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、フラッシュメモリ（ｆｌａｓｈ ｍｅｍｏｒｙ）等の半導体メモリ素子や、ハードディスク、光ディスク等によって実現される。記憶部１３０は、シンチレータ番号（Ｐ）、エネルギー値（Ｅ）、及び検出時間（Ｔ）が対応づけられた情報である計数情報、同時計数情報の通し番号であるコインシデンスＮｏ．に計数情報の組が対応づけられた同時計数情報、または同時計数情報を集計して得られる射影データ、再構成されたＰＥＴ画像等を記憶する。

処理回路１５０は、計測機能１５０ａ、特定機能１５０ｂ、再構成機能１５０ｃ、制御機能１５０ｄ、受付機能１５０ｅ、画像生成機能１５０ｆ、表示制御機能１５０ｇを有する。なお、計測機能１５０ａ及び特定機能１５０ｂの各機能については、後ほど詳しく説明する。

実施形態では、計測機能１５０ａ、特定機能１５０ｂ、再構成機能１５０ｃ、制御機能１５０ｄ、受付機能１５０ｅ、画像生成機能１５０ｆ、表示制御機能１５０ｇにて行われる各処理機能は、コンピュータによって実行可能なプログラムの形態で記憶部１３０へ記憶されている。処理回路１５０はプログラムを記憶部１３０から読み出し、実行することで各プログラムに対応する機能を実現するプロセッサである。換言すると、各プログラムを読み出した状態の処理回路１５０は、図１の処理回路１５０内に示された各機能を有することになる。

なお、図１においては単一の処理回路１５０にて、計測機能１５０ａ、特定機能１５０ｂ、再構成機能１５０ｃ、制御機能１５０ｄ、受付機能１５０ｅ、画像生成機能１５０ｆ、表示制御機能１５０ｇにて行われる処理機能が実現されるものとして説明するが、複数の独立したプロセッサを組み合わせて処理回路１５０を構成し、各プロセッサがプログラムを実行することにより機能を実現するものとしても構わない。換言すると、上述のそれぞれの機能がプログラムとして構成され、１つの処理回路１５０が各プログラムを実行する場合であってもよい。別の例として、特定の機能が専用の独立したプログラム実行回路に実装される場合であってもよい。

上記説明において用いた「プロセッサ」という文言は、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、ＧＰＵ（Ｇｒａｐｈｉｃａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）或いは、特定用途向け集積回路（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ：ＡＳＩＣ）、プログラマブル論理デバイス（例えば、単純プログラマブル論理デバイス（Ｓｉｍｐｌｅ Ｐｒｏｇｒａｍｍａｂｌｅ Ｌｏｇｉｃ Ｄｅｖｉｃｅ：ＳＰＬＤ）、複合プログラマブル論理デバイス（Ｃｏｍｐｌｅｘ Ｐｒｏｇｒａｍｍａｂｌｅ Ｌｏｇｉｃ Ｄｅｖｉｃｅ：ＣＰＬＤ）、及びフィールドプログラマブルゲートアレイ（Ｆｉｅｌｄ Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ：ＦＰＧＡ））等の回路を意味する。プロセッサは記憶部１３０に保存されたプログラムを読み出し実行することで機能を実現する。

なお、図１において、計測機能１５０ａ、特定機能１５０ｂ、再構成機能１５０ｃ、制御機能１５０ｄ、受付機能１５０ｅ、画像生成機能１５０ｆ、表示制御機能１５０ｇは、それぞれ、計測部、特定部、再構成部、制御部、受付部、生成部、表示制御部の一例である。

なお、処理回路１５０に代えて、フロントエンド回路１０２が、計測部、特定部等の処理を担ってもよい。

処理回路１５０は、再構成機能１５０ｃにより、フロントエンド回路１０２から取得したデータに基づいて、ＰＥＴ画像の再構成を行い、画像生成機能１５０ｆにより、画像の生成を行う。

処理回路１５０は、制御機能１５０ｄにより、架台装置１及びコンソール装置２を制御することによって、ＰＥＴ装置１００の全体制御を行う。例えば、処理回路１５０は、制御機能１５０ｄにより、ＰＥＴ装置１００における撮影を制御する。また、処理回路１０５は、制御機能１５０ｄにより、寝台駆動部１０６を制御する。

処理回路１５０は、受付機能１５０ｅにより、入力装置１１０を通じて、ユーザから情報の入力を受け付ける。また、処理回路１５０は、表示制御機能１５０ｇにより、ＰＥＴ画像その他のデータをディスプレイ１２０に表示させる。また、処理回路１５０は、画像生成機能１５０ｆにより、種々の画像の生成を行う。

続いて、実施形態に係る背景について簡単に説明する。

ＰＥＴ画像の再構成を行う際において、ＴｏＦ（Ｔｉｍｅ Ｏｆ Ｆｌｉｇｈｔ）カーネルを用いて画像再構成が行われる。ＴｏＦカーネルとは、２つの検出器で検出された信号の検出時間差の関数として表された、検出イベントの確率密度関数である。ここで、画像再構成に用いられるＴｏＦカーネルの時間幅が、検出器のＴｏＦ性能に整合するものでない場合、再構成後の画像の画質が低下する場合がある。例えば、検出器のＴｏＦ性能と比較して、ＴｏＦカーネルの時間幅が小さい場合には、画像再構成後の放射性物質の位置の推定位置が、本来の分布よりも局所的に集積してしまう。一方、検出器のＴｏＦ性能と比較して、ＴｏＦカーネルの時間幅が大きい場合には、画像再構成後の放射線物質の位置の推定位置が、本来の分布より周辺部分に分布してしまう。従って、検出器のＴｏＦ性能と、画像再構成に用いるＴｏＦカーネルとが、できるだけ一致することが望ましい。

しかしながら、例えば、イベント数が少ない場合等、観測揺らぎが大きくなり、ＴｏＦカーネルの算出に誤差が生じる場合がある。そこで、実施形態では、発光モデルに基づいて、同時計数イベント毎に確率分布モデルを生成したのち、それらの確率分布を加算することによりＴｏＦカーネルの算出を行う。このことで、イベントごとの観測値の平均からの揺らぎを考慮したＴｏＦカーネルを生成して画像再構成をすることができるので、例えば少ない同時計数イベントで観測揺らぎが大きく見えてしまう場合等でも、再構成後のＰＥＴ画像の画質を向上させることができる。

かかる背景に基づいて、実施形態に係る核医学診断装置は、処理回路１５０を備える。処理回路１５０は、計測機能１５０ａにより、第１の検出器で検出される第１のイベントに関する第１の光子数情報と、第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを計測する。処理回路１５０は、特定機能１５０ｂにより、第１の光子数情報と第２の光子数情報とに基づいて、第１のイベントと第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する。

また、実施形態に係るデータ処理方法は、第１の検出器で検出される第１のイベントに関する第１の光子数情報と、第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを計測し、第１の光子数情報と第２の光子数情報とに基づいて、第１のイベントと第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する。

また、実施形態に係るプログラムは、第１の検出器で検出される第１のイベントに関する第１の光子数情報と、前記第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを計測し、第１の光子数情報と第２の光子数情報とに基づいて、第１のイベントと第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する処理をコンピューターに実行させる。

以下、図２～図１０を用いて、実施形態に係る核医学診断装置が行う処理について、詳しく説明する。

図２に、実施形態に係る核医学診断装置としてのＰＥＴ装置１００が行う処理の流れが示されている。以下のフローチャートでは、簡単のため、第１の検出器ｌと第２の検出器ｍに係る確率分布モデル（ＴｏＦカーネル）を生成する場合についてのみ記載している。しかしながら、典型的には、画像再構成においては、すべてのｌとｍの組について処理回路１５０は確率分布モデルを生成した上でステップＳ７０の画像再構成処理を行う。この場合、処理回路１５０は、すべてのｌとｍの組についてステップＳ１０からステップＳ６０の処理を繰り返した上で、ステップＳ７０の処理を行う。

ステップＳ１０において、処理回路１５０は、処理を行う検出器の組を示すｌとｍの値を選択する。すなわち、処理回路１５０は、第１のイベントを検出する検出器である第１の検出器ｌ（ｌ番目の検出器）と、第２のイベントを検出する検出器である第２の検出器ｍ（ｍ番目の検出器）とを特徴づけるｌとｍの値を選択する。なお、ここでいう第１のイベントと第２のイベントとで、一つの同時計数イベントが構成される。別の言い方をすれば、同時計数イベントの各々は、第１のイベントと第２のイベントの２つのイベントからなる。

続いて、ステップＳ２０において、処理回路１５０は、同時計数イベントの番号を示すｉの値に１をセットする。ここで、ｉ=1とは、処理回路１５０が、１番目の同時計数イベントの処理を行うことを意味する。

続いて、ステップＳ３０において、処理回路１５０は、特定機能１５０ｂにより、ｌ番目の検出器で検出されたイベントと、ｍ番目の検出器で検出されたイベントとの同時計数イベントのうちのｉ番目の同時計数イベントに対して、コインシデンスに関する確率分布モデルｐ_{ｃｔ；ｌ、ｍ} ^{ｉ、（ｄｏｉ）}（θ_ｌ ^ｉ,θ_ｍ ^ｉ、ｔ）を算出する。ここで、θ_ｌ ^ｉ、θ_ｍ ^ｉは、ｉ番目の同時計数イベントにおいて、それぞれｌ番目の検出器及びｍ番目の検出器について取得されるパラメータを示す。また、（ｄｏｉ）の記号は、シンチレータ内での発光位置の不定性を考慮してもよいし、しなくてもよいことを意味し、シンチレータ内での発光位置の不定性を考慮する実施形態においては、処理回路１５０は、特定機能１５０ｂにより、シンチレータ内での発光位置の不定性を考慮したコインシデンスに関する確率分布モデルｐ_{ｃｔ；ｌ、ｍ} ^{ｉ、ｄｏｉ}（θ_ｌ ^ｉ,θ_ｍ ^ｉ、ｔ）を算出し、シンチレータ内での発光位置の不定性を考慮しない実施形態においては、シンチレータ内での発光位置の不定性を考慮しないコインシデンスに関する確率分布モデルｐ_{ｃｔ；ｌ、ｍ} ^ｉ（θ_ｌ ^ｉ,θ_ｍ ^ｉ、ｔ）を算出することを意味する。

図３に、図２のステップＳ３０の処理について、より詳細な処理の流れが示されている。すなわち、図３のステップＳ１００～ステップＳ５００は、図２のステップＳ３０をより詳細に説明したものとなる。

はじめに、ステップＳ１００において、処理回路１５０は、計測機能１５０ａにより、ｉ番目の同時計数イベントについて、第１の検出器ｌで検出される第１のイベントに関する第１の光子数情報を取得する。

ここで、第１の光子数情報とは、例えば、当該イベントについてのシンチレーション光子数Ｎ^ｉ _Ｓ;lである。また、別の例として、第１の光子数情報とは、例えば、当該イベントについてのチェレンコフ光子数Ｎ^ｉ _Ｃ;lである。また、別の例として、第１の光子数情報とは、当該イベントについてのシンチレーション光子数Ｎ^ｉ _Ｓ;l及びシンチレーション光子数Ｎ^ｉ _Ｓ;lである。すなわち、処理回路１５０は、計測機能１５０ａにより、シンチレーション光子数Ｎ^ｉ _Ｓ;l、チェレンコフ光子数Ｎ^ｉ _Ｃ;l、またはシンチレーション光子数Ｎ^ｉ _Ｓ;l及びチェレンコフ光子数Ｎ^ｉ _Ｃ;lを、第１の光子数情報として計測する。ここで、処理回路１５０が計測機能１５０ａにより取得する第１の光子数情報は、イベントごとに取得されるものである。

ここで、ＰＥＴ装置１００は、第１の検出器ｌで取得されるエネルギースペクトルに基づいて、第１の光子数情報を取得してもよい。

図４に、検出器３内での発光モデルと、処理回路１５０が計測機能１５０ａにより取得する光子スペクトルに関する説明が示されている。検出器３のシンチレータ１１にガンマ線１０が入射すると、ガンマ線１０はシンチレータ１１にエネルギーを光電吸収またはコンプトン散乱の形で付与するとともに、高エネルギー電子１２が発生する。曲線１６は、高エネルギー電子１２のエネルギーＥ_ｅを時間ｔの関数として示したものである。高エネルギー電子１２は、シンチレータ１１中を走行するとともにエネルギーを失っていくが、エネルギーが閾値Ｅ_ｔｈより大きい間は、チェレンコフ光１３を周囲に放出する。

曲線１７は、チェレンコフ光の光子密度ｎ_ｃ（ｔ）を、時間ｔのガウス関数で近似したものである。具体的には、チェレンコフ光の光子密度ｎ_ｃ（ｔ）は、Ｎ_ｃをチェレンコフ光の光子数、σ_Ｃをチェレンコフ光スペクトルの半値幅、ｔ_０をチェレンコフ光のピーク時刻として、以下の式（１）で与えられる。

一方、高エネルギー電子１２は、走行に伴い、シンチレータ１１の中で価電子１４を励起し、励起された価電子１４が基底状態に戻るとき、シンチレーション光１５が放出される。曲線１８は、価電子の励起密度ｎ_ｅｘ（ｔ）を、時刻ｔのガウス関数で近似したものである。価電子の励起密度ｎ_ｅｘ（ｔ）は、Ｎ_ｓをシンチレーション光の光子数、σ_Ｓをシンチレーション光スペクトルの半値幅、ｔ_１をシンチレーション光のピーク時刻として、以下の式（２）で与えられる。

曲線２１は、時刻２０において励起された価電子がもとの状態に遷移していく様子を示している。励起された価電子がもとの状態に遷移するのに伴ってシンチレーション光が発生する。曲線２２は、シンチレーション光の光子密度ｎ_ｓ（ｔ）を、時刻ｔの関数として示したものである。シンチレーション光の光子密度ｎ_ｓ（ｔ）は、緩和定数をτ_Sとして、価電子の励起密度ｎ_ｅｘ（ｔ）を用いて、以下の式（３）で与えられる。

ここで、価電子の励起密度ｎ_ｅｘ（ｔ）の関数形を式（２）で仮定した場合、式（３）に式（２）を代入すると、以下の式（４）が得られる。

すなわち、ＰＥＴ装置１００は、第１の検出器ｌ及び第２の検出器ｍの各々で取得されるエネルギースペクトルを取得し、当該エネルギースペクトルを基に曲線１７や曲線２２に示されるスペクトルを生成し、生成したスペクトルに基づいて、第１の光子数情報及び第２の光子数情報を取得してもよい。

なお、発光の分布関数ｐ_ｐｈ（ｔ）は、チェレンコフ光の光子密度ｎ_ｃ（ｔ）とシンチレーション光の光子密度ｎ_ｓ（ｔ）との和に比例するから、以下の式（５）が成り立つ。

ここで、発光の分布関数ｐ_ｐｈ（ｔ）は、以下の式（６）のように規格化されている。

なお、ステップＳ１００において、処理回路１５０は、計測機能１５０ａにより、光子数情報以外の情報を、第１の検出器ｌに係るパラメータセットとして、取得してもよい。一例として、処理回路１５０は、計測機能１５０ａにより、チェレンコフ光に対する光センサの時間分解能σ_Ｃを、第１の光子数情報に代えて、または第１の光子数情報に加えて、第１の検出器ｌにおいて取得するパラメータセットθ_ｌの一つとして、取得してもよい。また、別の例として、処理回路１５０は、計測機能１５０ａにより、シンチレーション光に対する光センサの時間分解能σ_Ｓを、第１の光子数情報に代えて、または第１の光子数情報に加えて、第１の検出器ｌにおいて取得するパラメータセットθ_ｌの一つとして、取得してもよい。また、別の例として、処理回路１５０は、計測機能１５０ａにより、過剰ノイズ発生確率を、第１の検出器ｌにおいて取得するパラメータセットθ_ｌの一つとして、取得してもよい。ここでいう過剰ノイズの発生確率とは、例えばオプティカルクロストークやアフターパルス等のＡＰＤ（Ａｖａｌａｎｃｈｅ ｐｈｏｔｏ ｄｉｏｄｅ）過剰ノイズの発生確率を言う。

ステップＳ２００において、処理回路１５０は、計測機能１５０ａにより、ｉ番目の同時計数イベントについて、第１の検出器ｌとは異なる第２の検出器ｍで検出される第２のイベントに関する第２の光子数情報を取得する。

ここで、第２の光子数情報とは、例えば、当該同時計数イベントについてのシンチレーション光子数Ｎ^ｉ _Ｓ;ｍである。また、別の例として、第２の光子数情報とは、例えば、当該同時計数イベントについてのチェレンコフ光子数Ｎ^ｉ _Ｃ;ｍである。また、別の例として、第２の光子数情報とは、当該同時計数イベントについてのシンチレーション光子数Ｎ^ｉ _Ｓ;ｍ及びシンチレーション光子数Ｎ^ｉ _Ｓ;ｍである。すなわち、処理回路１５０は、計測機能１５０ａにより、シンチレーション光子数Ｎ^ｉ _Ｓ;ｍ、チェレンコフ光子数Ｎ^ｉ _Ｃ;ｍ、またはシンチレーション光子数Ｎ^ｉ _Ｓ;ｍ及びチェレンコフ光子数Ｎ^ｉ _Ｃ;ｍを、第２の光子数情報として計測する。ここで、処理回路１５０が計測機能１５０ａにより計測する第２の光子数情報は、イベントごとに取得されるものである。

ＰＥＴ装置１００は、第２の検出器ｍで取得されるエネルギースペクトルに基づいて、第２の光子数情報を取得してもよい。

また、処理回路１５０は、取得機能１５０ａにより、チェレンコフ光に対する光センサの時間分解能σ_Ｃ、シンチレーション光に対する光センサの時間分解能σ_Ｓ、過剰ノイズ発生確率等を、第２の検出器ｍにおいて取得するパラメータセットθ_ｍとして、取得してもよい。

ステップＳ３００において、処理回路１５０は、特定機能１５０ｂにより、ステップＳ１００で取得された第１の光子数情報に基づいて、ｉ番目の同時計数イベントにおける、第１の検出器ｌに対する第１の発光確率モデルp_ｔ；ｌ ^{ｉ,(ｄｏｉ)}を算出する。ここで、処理回路１５０は、特定機能１５０ｂにより、第１の検出器ｌに対する第１の発光確率モデルp_ｔ；ｌ ^{ｉ,(ｄｏｉ)}を、同時計数イベントごとに定める。

図５にステップＳ３００の処理の一例がより詳細に説明されている。図５のフローチャートは、ステップＳ３００の処理の一例をより詳細に説明したフローチャートである。

ステップＳ３１０において、処理回路１５０は、特定機能１５０ｂにより、ステップＳ１００で取得された第１の検出器ｌにおける第１の光子数情報に基づいて、発光の分布関数モデルpⁱ _ｐｈ，ｌ (t)を生成する。発光の分布関数p_ｐｈ (t)は、時刻ｔに発光が起こる確率を示す分布関数である。発光の分布関数モデルpⁱ _ｐｈ,l (t)は、時刻に関して全時間で積分すると１になる。また、発光の分布関数モデルpⁱ _ｐｈ，ｌ(t)は、検出器ごと及び同時計数イベントごとに定められる。

処理回路１５０は、特定機能１５０ｂにより、第１の光子数情報に基づいて、例えば式（１）、（４）、（５）、（６）等を用いて発光の分布関数モデルpⁱ _ｐｈ,l (t)を生成する。第１の光子数情報が、検出器ごと及び同時計数イベントごとに取得されているので、処理回路１５０は、特定機能１５０ｂにより、発光の分布関数モデルp_ｐｈ,l ^ｉ (t)を、検出器ごと及び同時計数イベントごとに算出する。図６に、発光の分布関数p_ｐｈ(t)の一例が示されている。図６において、曲線３０は、発光の分布関数p_ｐｈ(t)を示している。また、検出時刻を定義するためのフォトンの閾値（Ｎ’）を併記する。

続いて、ステップＳ３２０Ａにおいて、処理回路１５０は、特定機能１５０ｂにより、発光の分布関数モデルpⁱ _ｐｈ,l (θⁱ _ｌ,t)に基づいて、第１の発光確率モデルｐⁱ _ｔ,l(θⁱ _ｌ, t)を算出する。

ここで、第１の発光確率モデルｐⁱ _ｔ,l(≧Ｎ’，θⁱ _ｌ, t)は、検出事象に含まれる複数の光子の内の最初のＮ’個以上の光子を検出する確率密度に関するモデルである。換言すると、第１の発光確率モデルｐ_ｔ,l(≧Ｎ’，t)は、検出事象に含まれる複数の光子の内の所定の数の光子を検出する確率密度に関するモデルである。図７に、第１の発光確率モデルｐ_ｔ,l(≧Ｎ’，t)が示されている。ここで、光子の検出閾値Ｎ’はある値に固定しているため、図示していない。図７のグラフ３１は、第１の発光確率モデルｐ_ｔ,l(t)の概形を示している。なお、第１の発光確率モデルｐ_ｔ,l(t)は、検出器ごと及び同時計数イベントごとに定められる。以下、必要に応じて、検出器の添え字、同時計数イベントの番号を表す添え字、及び取得されるパラメータの集合を表す引数θ^ｉ _ｌ等については適宜省略する。

検出事象に含まれる複数の光子の内の最初のＮ'個の光子を検出する確率密度である第１の発光確率モデルｐ_t（≧Ｎ’，ｔ）は、以下のように算出することができる。

はじめに、時刻ｔまでにＮ’個未満のフォトンしか検出しない確率Ｐ’（≧Ｎ’，ｔ）に関する方程式を導出する。これは、時刻ｔまでに高々Ｎ’－１個のフォトンしか検出していない確率であり、時刻ｔまでにｋ個のフォトンを検出している確率Ｐｄ（ｋ，ｔ）を用いて以下の式（７）で表される。

時刻ｔまでにｋ個のフォトンを検出している確率は、フォトンの検出事象を２項分布で評価することで、以下の式（８）のように微分形式で表現できる。

この微分方程式を順次ｋ＝０、１，２，…に対して解くと、帰納的に以下の式（９）が成り立つ。

式（７）より、時刻ｔまでにＮ’個未満のフォトンしか検出しない確率Ｐ’（≧Ｎ’，ｔ）は以下の式（１０）のように書くことが出来る。

これを式（８）に代入することで、検出事象に含まれる複数の光子の内の最初のＮ’個の光子を検出する確率密度である第１の発光確率モデルｐ_t（≧Ｎ’，ｔ）は、以下の式（１１）のように表すことができる。

ここで、関数F（ｔ）は発光の分布関数ｐ_ｐｈ（ｔ）の累積分布関数であり、以下の式（１２）のように書くことが出来る。

以下の実施例では簡単のため、初めてフォトンを検出した時刻で評価した確率分布、すなわち、Ｎ’＝０とした場合を例にとる。

すなわち、初の光子を検出する確率密度である第１の発光確率モデルをp_ｔ (t)とし、以下の式（１３）のように陽に表す。

すなわち、処理回路１５０は、特定機能１５０ｂにより、式（１２）及び式（１３）を用いることで、発光の分布関数ｐ_ph(t)に基づいて、第１の発光確率モデルp_ｔ(t)を計算することができる。このように、処理回路１５０は、特定機能１５０ｂにより、式（１２）及び式（１３）に第１の検出器における発光の分布関数モデルp_ｐｈ (t)を代入し、第１の光子数情報を式（１３）の右辺のＮに代入して式（１３）の右辺を評価することで、第１の光子数情報及び同時計数イベントごとに定められた第１の検出器における発光の分布関数モデルp_ｐｈ (t)に基づいて、第１のイベントについての第１の発光確率モデルｐ_ｔ(t)を算出する。

なお、上述の発光の分布関数ｐ_ph(t)及び第１の第１の発光確率モデルp_t(t)は、各検出器ごと及び同時計数イベントごとに定めることができる。従って、検出器の添え字及び同時計数イベントに関する添え字を明示的に記載すると、式（１３）は、以下の式（１４）のように書くこともできる。

ここで、p^ｉ _ｔ；ｌ(t)はｉ番目の同時計数イベントの検出器ｌにおける第１の発光確率モデル、ｐ^ｉ _ｐｈ；ｌ(t)はｉ番目の同時計数イベントの検出器ｌにおける発光の分布関数、Ｎ^ｉ _ｓ；ｌはｉ番目の同時計数イベントの検出器ｌにおけるシンチレーション光の光子数、Ｎ^ｉ _ｃ；ｌはｉ番目の同時計数イベントの検出器ｌにおけるチェレンコフ光の光子数を示す。

なお、第１の発光確率モデルを算出の際、処理回路１５０は、特定機能１５０ｂにより、シンチレータ内での発光位置の不定性を考慮した計算を行うこともできる。例えば、図８に示されるように、ガンマ線４０がシンチレータ１１の深さｘにある点４２で発光し、矢印４３で示されるように、シンチレーション光が物質中の光速ｃ／ｎで屈折率ｎの媒体中を伝播してシンチレータ１１の端４４に到達した場合、シンチレータ内の発光位置によって、ガンマ線４０がシンチレータに入射してからシンチレータ１１の端に到達するまでの時刻が異なるので、信号を検出する時刻が異なる。

すなわち、ガンマ線４０がシンチレータ１１に入射してから深さｘの地点まで到達するのに要する時刻ｔ_１は、以下のように、式（１５）で与えられる。

また、点４２で発光した光が深さｘからシンチレータ１１の端まで到達するのに要する時刻ｔ_２は、以下のように、式（１６）で与えられる。

従って、ガンマ線４０がシンチレータ１１に入射してからシンチレーション光がシンチレータ１１の端まで到達するのに要する時刻ｔ’は、以下のように、式（１７）で与えられる。

式（１７）より、シンチレータ内の発光位置ｘによって、信号を検出する時刻が異なることがわかる。

以上をもとに、式（１３）を基に、シンチレータ内の発光位置ｘを考慮した定式化を考えると、シンチレータ内の発光位置ｘを考慮した取り扱いを考える場合、以下の式（１８）が成り立つ。

ここで、p_ｔ、ｌ ^ｉ,dｏｉ(t)は、同時計数イベントごとに定められた、シンチレータ内での発光位置の不定性を考慮した第１の発光確率モデルであり、ｐ_ｐｈ、ｌ ^{ｉ,ｄｏｉ}(t)は、同時計数イベントごとに定められた、シンチレータの位置の発光位置の不定性を考慮した発光の分布関数モデルである。

なお、式（１８）は、処理回路１５０が計測機能１５０ａによりステップＳ１００において取得するパラメータθ_ｌを関数の引数として明示的に表現すると、以下の式（１９）及び式（２０）のようにも書くことができる。

ここで、θ^ｉ _ｌは、同時計数イベントごとに定められた、検出器ｌにおいて取得されるパラメータセットを示しており、例えば、同時計数イベントごとに定められたシンチレーション光子数Ｎ^ｉ _ｓ；ｌ、チェレンコフ光子数Ｎ^ｉ _ｃ；ｌ、チェレンコフ光に対する光センサの時間分解能σ^ｉ _Ｃ、シンチレーション光に対する光センサの時間分解能σ^ｉ _Ｓ等のパラメータの組を表す。

式（１８）に戻り、シンチレータの位置の発光位置の不定性を考慮した場合のステップＳ３００の処理について説明する。図９は、シンチレータ内での発光位置の不定性を考慮した場合の、発光確率モデルの算出の流れを示したフローチャートである。

図５で示した場合と同様に、ステップＳ３１０において、処理回路１５０は、特定機能１５０ｂにより、第１の検出器ｌにおける第１の光子数情報等のパラメータθ_ｌに基づいて、同時計数イベントごとの発光の分布関数モデルｐ_ph,l(t)を生成する。

続いて、ステップＳ３１１において、処理回路１５０は、特定機能１５０ｂにより、シンチレータ内での位置ｘごとに、発光確率密度ｐ_ｄｏｉ（ｘ）を、減弱係数μの値に基づいて算出する。位置ｘにおける発光確率密度ｐ_ｄｏｉ（ｘ）は、以下の式（２１）で表される。

なお、発光確率密度ｐ_ｄｏｉ（ｘ）は、以下の式（２２）に示されるように、規格化されている。

なお、図８において、グラフ４５には、減弱率μexp(-μx)が示されている。

図９に戻り、ステップＳ３１２において、処理回路１５０は、特定機能１５０ｂにより、ステップＳ３１０で生成された、同時計数イベントごとに定められた発光の分布関数である分布関数モデルｐ^ｉ _ｐｈ、ｌ（ｔ）と発光確率密度ｐ_ｄｏｉ（ｘ）とに基づいて、同時計数イベントごとに定められた発光位置不確定性を考慮した発光の分布関数モデルｐ_ｐｈ；ｌ ^i,doi(ｔ)を第１の発光確率モデルとして算出する。

具体的には、処理回路１５０は、特定機能１５０ｂにより、発光位置不確定性を考慮した発光の分布関数モデルｐ^{ｉ，ｄｏｉ} _ｐｈ；ｌ（ｔ）を、ステップＳ３１０で生成された分布関数モデルｐⁱ _ｐｈ、ｌ（ｔ）と、ステップＳ３１１で算出された発光確率密度ｐ_ｄｏｉ（ｘ）とに基づいて、次の式（２３）に基づいて算出する。

ここで、Ｌはシンチレータの長さであり、ｘは発光位置であり、ｔ’は式（１７）において説明した、ガンマ線がシンチレータに入射してからシンチレーション光がシンチレータの端まで到達するまでの遅延時間である。

続いて、ステップＳ３２０Ｂにおいて、処理回路１５０は、特定機能１５０ｂにより、同時計数イベントごとに定められた、発光位置不確定性を考慮した発光の分布関数モデルｐ^{i，ｄｏｉ} _ｐｈ；ｌ（ｔ）に基づいて、同時計数イベントごとに定められた、発光位置不確定性を考慮した発光確率モデルｐ^{ｉ,ｄｏｉ} _ｔ；ｌ（θ_l;t）を算出する。具体的には、式（２４）及び式（２５）に示されるように、処理回路１５０は、特定機能１５０ｂにより、式（１３）と同様に、同時計数イベントごとに定められた発光位置不確定性を考慮した発光の分布関数モデルｐ^i,doi _ｐｈ（ｔ）に基づいて、同時計数イベントごとに定められた発光位置不確定性を考慮した発光確率モデルｐ^i,dｏｉ _t（t）を算出する。

なお、式（２３）において、式（１７）を基に被積分関数のｘをｔ’に変数変換すると、以下の式（２６）が得られるので、処理回路１５０は、特定機能１５０ｂにより、式（２４）の右辺の評価において、式（２６）を用いて算出を行ってもよい。

このようにして、処理回路１５０は、特定機能１５０ｂにより、シンチレータ内での発光位置の不定性を考慮して、第１の発光確率モデルおよび第２の発光確率モデルを算出する。

図３に戻り、ステップＳ４００において、処理回路１５０は、特定機能１５０ｂにより、ステップＳ３００と同様の処理を行うことにより、ステップＳ２００で取得された第２の光子数情報に基づいて、第２の検出器ｍに対する、同時計数イベントごとに定められた第２の発光確率モデルp^{ｉ、（ｄｏｉ）} _ｔ；ｍを算出する。

具体的には、処理回路１５０は、特定機能１５０ｂにより、ステップＳ２００で取得された第２の検出器ｍにおける同時計数イベントごとに定められた第２の光子数情報に基づいて、例えば式（１）、（４）、（５）、（６）等に相当する式を用いて、同時計数イベントごとに定められた発光の分布関数モデルp^ｉ _ｐｈ；ｍ (t)を生成する。

続いて、処理回路１５０は、特定機能１５０ｂにより、第２の光子数情報及び第２の検出器ｍにおける発光の分布関数モデルp^ｉ _ｐｈ；ｍ(t)に基づいて、同時計数イベントごとに定められた第２の発光確率モデルｐ^ｉ _ｔ；ｍ(t)を、例えば式（１４）等に相当する式を用いて算出する。

また、第２の発光確率モデルを算出の際、処理回路１５０は、特定機能１５０ｂにより、シンチレータ内での発光位置の不定性を考慮した計算を、式（１７）及び式（２１）～（２６）と同様の式を用いて行うこともできる。この場合、処理回路１５０は、特定機能１５０ｂにより、シンチレータ内での位置ｘごとに、発光確率密度ｐ_ｄｏｉ（ｘ）を、減弱係数μの値に基づいて算出する。続いて、処理回路１５０は、特定機能１５０ｂにより、同時計数イベントごとに定められた分布関数モデルｐ^ｉ _ｐｈ；ｍ（ｔ）と発光確率密度ｐ_ｄｏｉ（ｘ）とに基づいて、同時計数イベントごとに定められた発光位置不確定性を考慮した発光の分布関数モデルｐ^{ｉ、ｄｏｉ} _ｐｈ；ｍ（ｔ）を算出する。続いて、処理回路１５０は、特定機能１５０ｂにより、同時計数イベントごとに定められた発光位置不確定性を考慮した発光の分布関数モデルｐ^{ｉ、ｄｏｉ} _ｐｈ；ｍ（ｔ）に基づいて、同時計数イベントごとに定められた発光位置不確定性を考慮した発光確率モデルｐ^{ｉ、ｄｏｉ} _ｔ；ｍ（t）を算出する。

このようにして、処理回路１５０は、特定機能１５０ｂにより、シンチレータ内での発光位置の不定性を考慮して、同時計数イベントごとに定められた第２の発光確率モデルｐ_ｔ,ｍ(t)を算出する。

続いて、図２に戻り、ステップＳ５００において、処理回路１５０は、特定機能１５０ｂにより、ステップS３００で算出された第１の発光確率モデルｐ_ｔ；ｌ ^ｄｏｉ（t_１）と算出された第２の発光確率モデルｐ_ｔ；ｍ ^ｄｏｉ（t_２）とに基づいて、同時計数イベントごとに特定された、第１のイベントと第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルｐ^{ｉ、（ｄｏｉ）} _ｃｔ;l,m(θ_l, θ_ｍ,t)を特定する。

図１０において、曲線５３は、第１の検出器ｌのシンチレータ５１における同時計数イベントごとに定められた第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _ｔ；ｌ（θ_l,t_l）を、遅延時間t_lの関数としてプロットしたものを示している。また、曲線５４は、第２の検出器ｍのシンチレータ５２における第２の発光確率モデルｐ^{ｉ，ｄｏｉ} _ｔ；ｍ（θ_ｍ,t_ｍ）を、遅延時間ｔ_ｍの関数としてプロットしたものを示している。処理回路１５０は、特定機能１５０ｂにより、検出器ｌにおける検出時刻ｔ_ｌにおける事象と、検出器ｍにおける検出時刻ｔ_ｍにおける積事象を、検出時間差の関数として表現した、同時計数イベントごとに定められた確率分布モデルｐ^{ｉ，ｄｏｉ} _ｃｔ;l,m(θ_l, θ_ｍ,t)を、特定する。曲線５６は、このような確率分布モデルｐ^{ｉ，ｄｏｉ} _ｃｔ;l,m(θ_l, θ_ｍ,t)の例を示している。

ここで、時間幅ｄｔの間における、検出器ｌにおける検出時刻ｔ_ｌにおける事象と、検出器ｍにおける検出時刻ｔ_ｍにおける積事象の数は、以下の式（２７）で与えられる。

ここで、θ^ｉ _ｌは、ｉ番目の同時計数イベントに係る、検出器ｌにおいて取得されるパラメータセットを示している。また、パラメータθを、より具体的な形で表現すると、以下の式（２８）が成り立つ。

ここで、Ｎ^ｉ _ｓ；ｌ、Ｎ^ｉ _ｓ；ｍは、ｉ番目の同時計数イベントに係る、検出器ｌ及びｍにおけるシンチレーション光子数であり、Ｎ^ｉ _ｃ；ｌ、Ｎ^ｉ _ｃ；ｍは、ｉ番目の同時計数イベントに係る、検出器ｌ及びｍにおけるチェレンコフ光子数を示す。

また、検出器ｌにおける検出時刻ｔ_ｌと、検出器ｍにおける検出時刻ｔ_ｍとの検出時間の差は、以下の式（２９）で与えられる。

従って、ｐ^{ｉ，ｄｏｉ} _ｃｔ;l,m(θ_l, θ_ｍ,t)は、以下の式（３０）で与えられる。

また、パラメータθを、より具体的な形で表現すると、以下の式（３１）が成り立つ。

以上のように、図２のステップＳ３０の処理により、処理回路１５０は、特定機能１５０ｂにより、第１の検出器ｌと第２の検出器ｍとによって規定されるＬＯＲの検出時間差に関する確率分布モデルｐ^{ｉ、ｄｏｉ} _ｃｔ;l,m(θ^ｉ _l, θ^ｉ _ｍ,t)を、ｉ番目の同時計数イベントについて算出する。処理回路１５０は、この処理を、全ての同時計数イベントについて算出し、それらを加算することで、最終的な、第１の検出器ｌと第２の検出器ｍとによって規定されるＬＯＲの検出時間差に関する確率分布モデルｐ^ｄｏｉ _ｃｔ;l,m(θ^ｉ _l, θ^ｉ _ｍ,t)を算出し、それを基に、データの再構成を行うことができる。

ステップＳ４０において、処理回路１５０は、特定機能１５０ｂにより、すべての同時計数イベントについて、ステップＳ３０の処理が終わったか否かを判定する。処理回路１５０が、特定機能１５０ｂにより、すべての同時計数イベントについて、ステップＳ３０の処理が終わっていないと判定した場合（ステップＳ４０ Ｎｏ）、処理はステップＳ５０に進み、ｉ＋１番目のデータについて、処理が行われる。一方、処理回路１５０が、特定機能１５０ｂにより、すべての同時計数イベントについて、ステップＳ３０の処理が終わったと判定した場合（ステップＳ４０ Ｙｅｓ）、処理はステップＳ６０に進む。

ステップＳ６０において、処理回路１５０は、特定機能１５０ｂにより、ステップＳ３０で算出された、各同時計数イベントごとの確率分布モデルｐ^{ｉ、ｄｏｉ} _ｃｔ;l,m(θ^ｉ _l, θ^ｉ _ｍ,t)を加算して、第１の検出器ｌと第２の検出器ｍにおけるコインシデンスに関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}を生成する。

具体的には、第１の検出器ｌと第２の検出器ｍにおけるコインシデンスに関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}は、以下の式（３２）で表される。

ここで、第１の検出器ｌと第２の検出器ｍにおけるコインシデンスに関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}は、後述のステップＳ７０において、ＴｏＦカーネルとして用いられることで、画像再構成に用いられる。すなわち、処理回路１５０は、特定機能１５０ｂにより、各イベントごとに特定された確率分布モデルｐ^{ｉ、ｄｏｉ} _ｃｔ;l,m(θ^ｉ _l, θ^ｉ _ｍ,t)を、各イベントについて加算することで、第１の検出器ｌと第２の検出器ｍにおけるコインシデンスに関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}を、第１の検出器ｌと第２の検出器ｍとの間のＴｏＦカーネルとして生成する。このように、処理回路１５０は、特定機能１５０ｂにより、第１の光子数情報と第２の光子数情報とに基づいて、第１のイベントと第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}を特定する。

ここで、実施形態では、処理回路１５０は、特定機能１５０ｂにより、ステップＳ３０において同時計数イベントごとに確率分布モデルｐ^{ｉ、ｄｏｉ} _ｃｔ;l,m(θ^ｉ _l, θ^ｉ _ｍ,t)を生成した上で、これらをステップＳ６０において加算して、最終的な確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}を生成している。このことにより、処理回路１５０は、特定機能１５０ｂにより、事象ごとの観測値の平均からの揺らぎを考慮したＴｏＦカーネルを生成することができる。この結果、処理回路１５０は、特定機能１５０ｂにより、例えば少ない事象で観測揺らぎが大きく見えてしまう場合等でも、妥当なＴｏＦカーネルを生成することができる。

続いて、ステップＳ７０において、処理回路１５０は、再構成機能１５０ｃにより、ステップＳ６０で特定された、第１の光子数情報と第２の光子数情報とに基づいて規定されたコインシデンスに関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}をＴＯＦカーネルとして利用することで、第１の検出器ｌと第２の検出器ｍとの間のＬＯＲに対応するデータの再構成を行う。

図１１～図１３に、実施形態に係る手法の検証結果が示されている。この検証では、リストモードデータから、本手法を用いて確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}を生成し、実測値と比較した。

図１１に、ＢＧＯ結晶を用いて測定されたリストモードの散布図が示されており、図中の丸一つ一つ、例えばイベント６０及びイベント６１が一つのイベントに対応する。処理回路１５０は、特定機能１５０ｂにより、イベント６０を基に、イベント６０に対応する、イベントごとの確率分布モデルｐ_ｃｔ（ｔ）を算出した。曲線６３は、イベント６０に対応するｐ_ｃｔ（ｔ）である。また、処理回路１５０は、特定機能１５０ｂにより、イベント６１を基に、イベント６１に対応する、イベントごとの確率分布モデルｐ_ｃｔ（ｔ）を算出した。曲線６２は、イベント６１に対応するｐ_ｃｔ（ｔ）である。処理回路１５０は、特定機能１５０ｂにより、同様な処理をすべてのイベントについて行い、イベントごとの確率分布モデルｐ_ｃｔ（ｔ）を加算することにより、確率分布モデルｐ^ｔｏｔ（ｔ）を生成した。

図１２は、ＬＹＳＯにおける実測値と確率分布モデルｐ^ｔｏｔ _ｃｔ（ｔ）との比較を示している。曲線６５が実測値であり、曲線６４が、計算された確率分布モデルｐ^ｔｏｔ _ｃｔ（ｔ）である。また、図１３は、ＢＧＯにおける実測値と確率分布モデルｐ^ｔｏｔ _ｃｔ（ｔ）との比較を示している。曲線６６が実測値であり、曲線６７が、計算された確率分布モデルｐ^ｔｏｔ _ｃｔ（ｔ）である。どちらの場合でも、実測値と確率分布モデルｐ^ｔｏｔ（ｔ）とが概ね一致しており、同時計数イベントごとに発光モデルをもとに確率分布モデルｐ_ｃｔ（ｔ）を生成し、それらを加算して確率分布モデルｐ^ｔｏｔ _ｃｔ（ｔ）を生成する実施形態の手法の妥当性が確認された。

実施形態は、上述の例に限られない。上述の場合においては、検出器３がモノリシック検出器でない通常の検出器である場合について説明したが、実施形態はこれらに限られない。実施形態は、検出器３が、単結晶として構成され、１つのシンチレータ結晶が区分けされずに一体化したまま使用されるモノリシックシンチレータである場合であってもよい。当該シンチレータは、例えば、ＬＹＳＯ（Lutetium Yttrium Oxyorthosilicate）、ＬＳＯ（Lutetium Oxyorthosilicate）、ＬＧＳＯ（Lutetium Gadolinium Oxyorthosilicate）等の、シンチレータ結晶によって形成される。

検出器３がモノリシック検出器である場合、検出器３は一つのシンチレータ結晶が区分けされずに一体化したまま使用されるものであるが、検出器３が複数の仮想的な検出器からなると考えて、それぞれの仮想的な検出器の間のＴｏＦカーネルを生成することにより、同様の処理を行うことができる。例えば、検出器３が複数の検出素子アレイからなる場合、それぞれの検出素子アレイを１つの検出器と考えて、図２のフローチャートで示された処理と同様の処理を行うことができる。一例として、図２のステップＳ１０からステップＳ７０や図３のステップＳ１００からステップＳ３００、式（２０）、（３０）等における自然数ｌ及びｍを、検出器に付された番号の代わりに、検出素子アレイに付された番号であると読み替えて同様な処理を行うことで、処理回路１５０は、特定機能当該検出素子アレイの組み合わせに応じた検出時間差に関する確率分布モデルを特定することができる。

すなわち、検出器３がモノリシック検出器である場合、実施形態における核医学診断装置は、処理回路１５０を備える。処理回路１５０は、計測機能１５０ａにより、少なくとも第１の検出素子アレイｌと第２の検出素子アレイｍとが設けられたモノリシック検出器と、第１の検出素子アレイｌにて検出される第１のイベントに対応する第１の光子数情報と、第２の検出素子アレイｍにて検出され第１のイベントに対応する第２の光子数情報とを計測する。処理回路１５０は、特定機能１５０ｂにより、第１の光子数情報と第２の光子数情報とに基づいて、第１の検出素子アレイｌと第２の検出素子アレイｍとの組み合わせに応じた検出時間差に関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}を特定する。

また、実施形態に係るデータ処理方法は、少なくとも第１の検出素子アレイｌと第２の検出素子アレイｍとが設けられたモノリシック検出器の第１の検出素子アレイｌにて検出される第１のイベントに対応する第１の光子数情報と、第２の検出素子アレイｍにて検出され第１のイベントに対応する第２の光子数情報とを計測し、第１の光子数情報と第２の光子数情報とに基づいて、第１の検出素子アレイｌと第２の検出素子アレイｍとの組み合わせに応じた検出時間差に関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}を特定する。

また、実施形態に係るプログラムは、少なくとも第１の検出素子アレイｌと第2の検出素子アレイｍとが設けられたモノリシック検出器の第１の検出素子アレイｌにて検出される第１のイベントに対応する第１の光子数情報と、第２の検出素子アレイｍにて検出され第１のイベントに対応する第２の光子数情報とを計測し、第１の光子数情報と第２の光子数情報とに基づいて、第１の検出素子アレイｌと第２の検出素子アレイｍとの組み合わせに応じた検出時間差に関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}を特定する処理をコンピューターに実行させる。

以上のように、実施形態に係る核医学診断装置においては、処理回路１５０は、シンチレータ内の発光モデルに基づいて、検出時間差に関する確率分布モデルを、同時計数イベントごとに算出した上でそれらを加算して検出時間差に関する確率分布モデルを生成した。これにより、事象ごとの観測値の平均からの揺らぎを考慮したＴｏＦカーネルを生成して画像再構成をすることができるので、例えば少ない事象で観測揺らぎが大きく見えてしまう場合等でも、再構成後のＰＥＴ画像の画質を向上させることができる。

以上の実施形態に関し、発明の一側面および選択的な特徴として以下の付記を開示する。
(付記１)
本発明の一つの側面において提供される核医学診断装置は、計測部と、特定部とを備える。計測部は、第１の検出器で検出される第１のイベントに関する第１の光子数情報と、前記第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを計測する。特定部は、前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１のイベントと前記第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する。

（付記２）
前記コインシデンスに関する確率分布モデルをＴｏＦ（Ｔｉｍｅ Ｏｆ Ｆｌｉｇｈｔ）カーネルとして利用することで、前記第１の検出器と前記第２の検出器との間のＬＯＲ（Ｌｉｎｅ Ｏｆ Ｒｅｓｐｏｎｓｅ）に対応するデータの再構成を行う再構成処理部を更に備えてもよい。

（付記３）
前記計測部は、前記第１の光子数情報及び前記第２の光子数情報として、シンチレーション光子数を計測してもよい。

（付記４）
前記計測部は、前記第１の光子数情報及び前記第２の光子数情報として、チェレンコフ光子数を更に計測してもよい。

（付記５）
前記特定部は、各イベントごとに特定された前記確率分布モデルを加算して、前記第１の検出器と前記第２の検出器との間のＴｏＦカーネルを生成してもよい。

（付記６）
前記特定部は、前記第１の光子数情報に基づいて、前記第１のイベントについての第１の発光確率モデルを生成し、前記第２の光子数情報に基づいて、前記第２のイベントについての第２の発光確率モデルを生成し、前記第１の発光確率モデルと前記第２の発光確率モデルとに基づいて、前記確率分布モデルを特定してもよい。

（付記７）
前記特定部は、イベントごとに定められた発光分布関数モデルに基づいて、前記第１の発光確率モデルおよび前記第２の発光確率モデルを生成してもよい。

（付記８）
前記計測部により計測される前記第１の光子数情報及び前記第２の光子数情報は、イベントごとに取得されるものであってもよい。

（付記９）
前記特定部は、シンチレーター内での発光位置の不定性を考慮して、前記第１の発光確率モデルおよび前記第２の発光確率モデルを算出してもよい。

（付記１０）
前記特定部は、前記シンチレーター内での位置ごとに発光確率密度を減弱係数に基づいて算出し、発光の分布関数及び前記発光確率密度に基づいて、前記第１の発光確率モデルおよび前記第２の発光確率モデルを算出してもよい。

（付記１１）
前記第１の発光確率モデルは、検出事象に含まれる複数の光子の内の所定の数の光子を検出する確率密度に関するモデルであってもよい。

（付記１２）
本発明の一つの側面において提供されるデータ処理方法は、第１の検出器で検出される第１のイベントに関する第１の光子数情報と、前記第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを計測し、前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１のイベントと前記第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する。

（付記１３）
本発明の一つの側面において提供されるプログラムは、第１の検出器で検出される第１のイベントに関する第１の光子数情報と、前記第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを計測し、前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１のイベントと前記第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する処理をコンピューターに実行させる。

（付記１４）
本発明の一つの側面において提供される核医学診断装置は、計測部と、特定部とを備える。計測部は、少なくとも第１の検出素子アレイと第２の検出素子アレイとが設けられたモノリシック検出器と、前記第１の検出素子アレイにて検出される第１のイベントに対応する第１の光子数情報と、第２の検出素子アレイにて検出され前記第１のイベントに対応する第２の光子数情報とを計測する。特定部は、前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１の検出素子アレイと前記第２の検出素子アレイとの組み合わせに応じた検出時間差に関する確率分布モデルを特定する。

（付記１５）
本発明の一つの側面において提供されるデータ処理方法は少なくとも第１の検出素子アレイと第２の検出素子アレイとが設けられたモノリシック検出器の前記第１の検出素子アレイにて検出される第１のイベントに対応する第１の光子数情報と、第２の検出素子アレイにて検出され前記第１のイベントに対応する第２の光子数情報とを計測し、前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１の検出素子アレイと前記第２の検出素子アレイとの組み合わせに応じた検出時間差に関する確率分布モデルを特定する。

（付記１６）
本発明の一つの側面において提供されるプログラムは、少なくとも第１の検出素子アレイと第２の検出素子アレイとが設けられたモノリシック検出器の前記第１の検出素子アレイにて検出される第１のイベントに対応する第１の光子数情報と、第２の検出素子アレイにて検出され前記第１のイベントに対応する第２の光子数情報とを計測し、前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１の検出素子アレイと前記第２の検出素子アレイとの組み合わせに応じた検出時間差に関する確率分布モデルを特定する処理をコンピューターに実行させる。

以上説明した少なくとも１つの実施形態によれば、画質を向上させることができる。

いくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更、実施形態同士の組み合わせを行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。

１５０ 処理回路
１５０ａ 計測機能
１５０ｂ 特定機能
１５０ｃ 再構成機能
１５０ｄ 制御機能
１５０ｅ 受付機能
１５０ｆ 画像生成機能
１５０ｇ 表示制御機能

Claims (15)

1. 第１の検出器で検出される第１のイベントに関する第１の光子数情報と、前記第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを取得する計測部と、
   前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１のイベントと前記第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する特定部と、
   を備える、核医学診断装置。
2. 前記コインシデンスに関する確率分布モデルをＴｏＦ（Ｔｉｍｅ Ｏｆ Ｆｌｉｇｈｔ）カーネルとして利用することで、前記第１の検出器と前記第２の検出器との間のＬＯＲ（Ｌｉｎｅ Ｏｆ Ｒｅｓｐｏｎｓｅ）に対応するデータの再構成を行う再構成処理部を更に備える、請求項１に記載の核医学診断装置。
3. 前記計測部は、前記第１の光子数情報及び前記第２の光子数情報として、シンチレーション光子数を計測する、請求項１に記載の核医学診断装置。
4. 前記計測部は、前記第１の光子数情報及び前記第２の光子数情報として、チェレンコフ光子数を更に計測する、請求項３に記載の核医学診断装置。
5. 前記特定部は、イベントごとに特定された前記確率分布モデルを加算して、前記第１の検出器と前記第２の検出器との間のＴｏＦカーネルを生成する、請求項１に記載の核医学診断装置。
6. 前記特定部は、前記第１の光子数情報に基づいて、前記第１のイベントについての第１の発光確率モデルを生成し、前記第２の光子数情報に基づいて、前記第２のイベントについての第２の発光確率モデルを生成し、前記第１の発光確率モデルと前記第２の発光確率モデルとに基づいて、前記確率分布モデルを特定する、請求項１に記載の核医学診断装置。
7. 前記特定部は、イベントごとに定められた発光分布関数モデルに基づいて、前記第１の発光確率モデルおよび前記第２の発光確率モデルを生成する、請求項６に記載の核医学診断装置。
8. 前記計測部により計測される前記第１の光子数情報及び前記第２の光子数情報は、イベントごとに取得されるものである、請求項１に記載の核医学診断装置。
9. 前記特定部は、シンチレータ内での発光位置の不定性を考慮して、前記第１の発光確率モデルおよび前記第２の発光確率モデルを算出する、請求項６に記載の核医学診断装置。
10. 前記特定部は、前記シンチレータ内での位置ごとに発光確率密度を減弱係数に基づいて算出し、発光の分布関数及び前記発光確率密度に基づいて、前記第１の発光確率モデルおよび前記第２の発光確率モデルを算出する、請求項９に記載の核医学診断装置。
11. 第１の検出器で検出される第１のイベントに関する第１の光子数情報と、前記第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを取得し、
    前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１のイベントと前記第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する、データ処理方法。
12. 第１の検出器で検出される第１のイベントに関する第１の光子数情報と、前記第１の検出器とは異なる第２の検出器で検出される第２のイベントに関する第２の光子数情報とを取得し、
    前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１のイベントと前記第２のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルを特定する処理をコンピューターに実行させる、プログラム。
13. 少なくとも第１の検出素子アレイと第２の検出素子アレイとが設けられたモノリシック検出器と、前記第１の検出素子アレイにて検出される第１のイベントに対応する第１の光子数情報と、第２の検出素子アレイにて検出され前記第１のイベントに対応する第２の光子数情報とを取得する計測部と、
    前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１の検出素子アレイと前記第２の検出素子アレイとの組み合わせに応じた検出時間差に関する確率分布モデルを特定する特定部と、
    を備える、核医学診断装置。
14. 少なくとも第１の検出素子アレイと第２の検出素子アレイとが設けられたモノリシック検出器の前記第１の検出素子アレイにて検出される第１のイベントに対応する第１の光子数情報と、第２の検出素子アレイにて検出され前記第１のイベントに対応する第２の光子数情報とを取得し、
    前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１の検出素子アレイと前記第２の検出素子アレイとの組み合わせに応じた検出時間差に関する確率分布モデルを特定する、データ処理方法。
15. 少なくとも第１の検出素子アレイと第２の検出素子アレイとが設けられたモノリシック検出器の前記第１の検出素子アレイにて検出される第１のイベントに対応する第１の光子数情報と、第２の検出素子アレイにて検出され前記第１のイベントに対応する第２の光子数情報とを取得し、
    前記第１の光子数情報と前記第２の光子数情報とに基づいて、前記第１の検出素子アレイと前記第２の検出素子アレイとの組み合わせに応じた検出時間差に関する確率分布モデルを特定する処理をコンピューターに実行させる、プログラム。

Images