JP2024052403A - 核医学診断装置、データ処理方法及びプログラム - 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として、価電子の励起密度ｎ_ｅｘ（ｔ）を用いて、以下の式（３）で与えられる。

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

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

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

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

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

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

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

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

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

ステップＳ１００に続き、ステップＳ３００において、処理回路１５０は、算出機能１５０ｂにより、ステップＳ１００で取得された第１の光子数情報に基づいて、ｉ番目の同時計数イベントにおける、第１の検出器ｌに対する第１の発光確率モデルp_ｔ；ｌ ^{ｉ,(ｄｏｉ)}を算出する。なお、処理回路１５０は、算出機能１５０ｂにより、第１の検出器ｌに対する第１の発光確率モデル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(Δｔ)が示されている。ここで、光子の検出閾値Ｎ’はある値に固定しているため、図示していない。図７の曲線３１は、第１の発光確率モデルｐ_ｔ,l(Δｔ)の概形を示している。なお、第１の発光確率モデルｐ_ｔ,l(Δｔ)は、検出器ごと及び同時計数イベントごとに定められる。以下、必要に応じて、検出器の添え字、同時計数イベントの番号を表す添え字、及び取得されるパラメータの集合を表す引数θ^ｉ _ｌ等については適宜省略する。なお、曲線３２は、第１の発光確率モデルｐ_ｔ,l(t)の累積分布Ｐ（ｔ）を表す曲線である。曲線３２については後述する。

続いて、第１の発光確率モデルｐ_ｔ,l(≧Ｎ’，t)の算出方法について説明する。検出事象に含まれる複数の光子の内の最初のＮ’個以上の光子を検出する確率密度である第１の発光確率モデルｐ_t（≧Ｎ’，ｔ）は、以下のように算出することができる。

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

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

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

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

これを式（８）に代入することで、検出事象に含まれる複数の光子の内の最初のＮ’個の光子を検出する確率密度である第１の発光確率モデルｐ_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)を算出する。

続いて、ステップＳ４１０において、処理回路１５０は、特定機能１５０ｃにより、第１の発光確率モデルに基づいて、発光の検出確率が所定の閾値以上となる第１のタイミングを特定する。一例として、処理回路１５０は、特定機能１５０ｃにより、ステップＳ３００において算出された第１の発光確率モデルに基づいて、処理回路１５０は、特定機能１５０ｃにより、イベントごとの推定検出遅延時間Δｔ^ｌ、ｉ _ｐを、第１のタイミングとして特定する。

ここで、推定検出遅延時間とは、例えばＸ線がシンチレータに入射してから、発光の信号が検出されるまでにかかる時刻の推定値である。ここで、処理回路１５０は、特定機能１５０ｃにより、検出器ごと、同時計数イベントごとに、推定検出遅延時間の推定を行うことができる。処理回路１５０は、推定検出遅延時間を推定して検出時刻から差し引くことで、ＴｏＦスペクトルの先鋭化を実現することができ、画質を向上させることができる。

第１のタイミングの特定について、図７に戻り説明する。発光の検出確率が所定の閾値以上となる第１のタイミングを特定する第１の方法として、処理回路１５０は、特定機能１５０ｃにより、第１の発光確率モデルｐ^{ｉ,ｄｏｉ} _ｔ；ｌ（θ_l;t）により示される分布の累積分布P（ｔ）に基づいて、発光の検出確率が所定の閾値以上となる第１のタイミングを特定する。ここで、第１の発光確率モデルｐ^{ｉ,ｄｏｉ} _ｔ；ｌ（θ_l;t）により示される分布の累積分布P（ｔ）は、例えば、以下の式（２７）で与えられる。

ここで、処理回路１５０は、特定機能１５０ｃにより、第１の発光確率モデルｐ^{ｉ,ｄｏｉ} _ｔ；ｌ（θ_l;t）により示される分布の累積分布P（ｔ）が、所定の閾値ｐに達するような時刻、すなわち曲線３２と直線３３との交点となる時刻３４を、第１のタイミングΔｔ_ｐ ^ｌ，ｉとして特定する。これを式で表現すると、第１のタイミングΔｔ_ｐ ^ｌ，ｉは、以下の式（２８）で与えられる。

また、発光の検出確率が所定の閾値以上となる第１のタイミングを特定する第２の方法として、処理回路１５０は、特定機能１５０ｃにより、第１の発光確率モデルｐ^{ｉ,ｄｏｉ} _ｔ；ｌ（θ_l;t）により示される分布が極値を取る時刻３５を、第１のタイミングΔｔ_ｐ ^ｌ，ｉとして特定する。これを式で表現すると、第１のタイミングΔｔ_ｐ ^ｌ，ｉは、以下の式（２９）で与えられる。

処理回路１５０は、特定機能１５０ｃにより、第１のタイミングΔｔ_ｐ ^ｌ，ｉを、典型的には、同時計数イベントごと及び検出器ごとに特定する。しかしながら、実施形態はこれに限られず、処理回路１５０は、特定機能１５０ｃにより、イベント単位以外で、または検出器単位以外で、発光の検出確率が所定の閾値以上となる第１のタイミングを特定してもよい。

図３に戻り、ステップＳ４１０と同様にして、ステップＳ４２０において、処理回路１５０は、特定機能１５０ｃにより、第２の発光確率モデルに基づいて、発光の検出確率が所定の閾値以上となる第２のタイミングを特定する。一例として、処理回路１５０は、特定機能１５０ｃにより、ステップＳ３００において算出された第２の発光確率モデルに基づいて、処理回路１５０は、特定機能１５０ｃにより、イベントごとの推定検出遅延時間Δｔ^ｍ、ｉ _ｐを、第２のタイミングとして特定する。

発光の検出確率が所定の閾値以上となる第２のタイミングを特定する第１の方法として、処理回路１５０は、特定機能１５０ｃにより、第２の発光確率モデルｐ^{ｉ,ｄｏｉ} _ｔ；ｍ（θ_ｍ;t）により示される分布の累積分布P（ｔ）に基づいて、発光の検出確率が所定の閾値以上となる第１のタイミングを特定する。一例として、処理回路１５０は、特定機能１５０ｃにより、第１の発光確率モデルｐ^{ｉ,ｄｏｉ} _ｔ；ｍ（θ_ｍ;t）により示される分布の累積分布P（ｔ）が、所定の閾値ｐに達するような時刻を、第２のタイミングΔｔ_ｐ ^ｍ，ｉとして特定する。

また、発光の検出確率が所定の閾値以上となる第２のタイミングを特定する第２の方法として、処理回路１５０は、特定機能１５０ｃにより、第２の発光確率モデルｐ^{ｉ,ｄｏｉ} _ｔ；ｍ（θ_ｍ;t）により示される分布が極値を取る時刻を、第２のタイミングΔｔ_ｐ ^ｍ，ｉとして特定する。

ステップＳ４５０において、処理回路１５０は、補正機能１５０ｅにより、ステップＳ４１０で推定された第１のタイミングΔｔ_ｐ ^ｌ，ｉに基づいて、検出器ｌで検出されたイベントであってｉ番目の同時計数イベントと対応付けられているイベントについて、検出時刻の補正を行う。ここで、補正後の検出時刻ｔ^ｌ，ｉ _corrected は、補正前の検出時刻ｔ^ｌ，ｉ _dと、ステップＳ４１０で特定された第１のタイミングΔｔ_ｐ ^ｌ，ｉとを用いて、以下の式（３０）で表される。

換言すると、処理回路１５０は、測定機能１５０ｄにより、第１の検出器ｌで検出されるイベントの検出タイミングｔ^ｌ，ｉ _dを測定し、補正機能１５０ｅに基づいて、第１のタイミングΔｔ_ｐ ^ｌ，ｉに基づいて、検出タイミングｔ^ｌ，ｉ _dを補正し、補正後の検出時刻ｔ^ｌ，ｉ _correctedを得る。

また、これに伴い、処理回路１５０は、補正機能１５０ｅにより、検出時刻が補正された第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _{ｔ；ｌ；corrected}(θ_ｍ;t）を算出する。ここで、補正された第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _{ｔ；ｌ；corrected}(θ_ｌ;t）と、補正前の第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _{ｔ；ｌ；corrected}(θ_ｌ;t）との関係は、例えば以下の式（３１）で与えられる。

なお、説明を簡単にするために、ステップＳ５００以降の処理において、補正された第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _{ｔ；ｌ；corrected}(θ_ｌ;t）を、単に、第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _ｔ；ｌ(θ_ｌ;t）と記載するが、以降のコインシデンスに関する確率分布モデルを算出する処理において、第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _ｔ；ｌ(θ_ｌ;t）は、補正された第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _{ｔ；ｌ；corrected}(θ_ｌ;t）を意味する。

ステップＳ４６０において、ステップＳ４５０と同様の処理を行い、処理回路１５０は、補正機能１５０ｅにより、ステップＳ４００で推定された第２のタイミングΔｔ_ｐ ^ｍ，ｉに基づいて、検出器ｍで検出されたイベントであってｉ番目の同時計数イベントと対応付けられているイベントについて、検出時刻の補正を行う。すなわち、処理回路１５０は、測定機能１５０ｄにより、第１の検出器ｌで検出されるイベントの検出タイミングｔ^ｍ，ｉ _dを測定し、補正機能１５０ｅに基づいて、第２のタイミングΔｔ_ｐ ^ｍ，ｉに基づいて、検出タイミングｔ^ｍ，ｉ _dを補正し、補正後の検出時刻ｔ^ｍ，ｉ _correctedを得る。

また、これに伴い、処理回路１５０は、補正機能１５０ｅにより、検出時刻が補正された第２の発光確率モデルｐ^{ｉ，ｄｏｉ} _{ｔ；ｍ；corrected}(θ_ｍ;t）を算出する。ステップＳ５００以降の処理において、補正された第２の発光確率モデルｐ^{ｉ，ｄｏｉ} _{ｔ；ｍ；corrected}(θ_ｍ;t）を、単に、第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _ｔ；ｍ(θ_ｍ;t）と記載するが、ステップＳ５００以降のコインシデンスに関する確率分布モデルを算出する処理において、第２の発光確率モデルｐ^{ｉ，ｄｏｉ} _ｔ；ｍ(θ_ｍ;t）は、補正された第２の発光確率モデルｐ^{ｉ，ｄｏｉ} _{ｔ；ｍ；corrected}(θ_ｍ;t）を意味する。

続いて、ステップＳ５００において、処理回路１５０は、算出機能１５０ｂにより、ステップＳ３００で算出された第１の発光確率モデルｐ_ｔ；ｌ ^ｄｏｉ（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)を加算して、第１の検出器ｌと第２の検出器ｍにおけるコインシデンスに関する確率分布モデルｐ^ｔｏｔ _{ｃｔ；ｌ、ｍ}を生成する。

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

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

図１１～図１４に、実施形態に係る手法の検証結果が示されている。この検証では、第１のタイミングΔｔ_ｐ ^ｌ，ｉを算出するための閾値ｐの値をさまざまに変化させ、再構成画像において放射線物質の分布のピークが先鋭化するかどうかを確認した。

図１１に、事象１、２、３の３つの事象について、第１の発光確率モデルｐ^{ｉ，ｄｏｉ} _ｔ；ｌ(θ_ｌ;t）が示されている。この３つの事象について、式（２８）の右辺の閾値ｐを変化させて、それぞれの事象についての第１のタイミングΔｔ_ｐを算出した。直線６０がｐ＝０．１に対応し、直線６１がｐ＝０．４に対応し、直線６２がｐ＝０．８に対応する。また、矢印６３は、ｐ＝０．４の時の第１の事象に対する第１のタイミングΔｔ_{p_１}を、矢印６４は、ｐ＝０．４の時の第２の事象に対する第１のタイミングΔｔ_{p_２}を、矢印６５は、ｐ＝０．４の時の第３の事象に対する第１のタイミングΔｔ_{p_３}を示している。

図１２に、ｐ＝０．４及びｐ＝０．１の場合における、補正後のコインシデンススペクトルが示されている。グラフ７１はｐ＝０．４の時の補正後のコインシデンススペクトルであり、グラフ７０はｐ＝０．１の時の補正後のコインシデンススペクトルを示す。また、遅延時間が０付近でのカウント数を比較すると、ｐ＝０．１においてはカウント数が多いが幅も大きく、これに対して、ｐ＝０．４においては、カウント数も多く幅も小さくなった。従って、閾値をｐ＝０．４付近にした場合、コインシデンススペクトルにおいてカウント数が大きくなりまた幅も小さくなる。

また、実施形態に係る方法による補正を行った場合と行わない場合とで、再構成画像において、放射線物質の分布のピークが先鋭化するか否かを検証した。まず、図１３に示されるように、実測データ８０に対して、フィティングカーブ８１を生成し、ＴｏＦカーネルのスペクトルの同定を行った。続いて、得られたＴｏＦカーネルのスペクトルを基に、ＭＬ－ＥＭ法で再構成を行った。図１４に、再構成結果が示されている。グラフ９０は、実施形態に係る方法による補正を行った場合の再構成結果、グラフ９１は、実施形態に係る方法をによる補正を行わなかった場合の再構成結果を示している。グラフ９０とグラフ９１とを比較すると、ピーク付近でのカウント数が約２倍になっており、放射線物質の分布のピークが先鋭化することが確認できた。

以上のように、実施形態では、発光モデルに基づいて、ガンマ線がシンチレータに入射してから発光が起こるまでの遅延時刻を推定し、推定結果に基づいてＴｏＦカーネルの補正を行う。このことで、ＴｏＦスペクトルを先鋭化することが可能となり、再構成後のＰＥＴ画像の画質を向上させることができる。

以上の実施形態に関し、発明の一側面および選択的な特徴として以下の付記を開示する。

(付記１)
本発明の一つの側面において提供される核医学診断装置は、取得部と、算出部と、特定部と、測定部と、補正部とを備える。取得部は、第１の検出器で検出される第１の光子数情報を取得する。算出部は、前記第１の光子数情報に基づいて前記第１の検出器に対応する第１の発光確率モデルを算出する。特定部は、前記第１の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第１のタイミングを特定する。測定部は、前記第１の検出器で検出されるイベントの検出タイミングを測定する。補正部は、前記第１のタイミングに基づいて、前記検出タイミングを補正する。

(付記２)
前記第１の光子数情報は、シンチレーション光子数を含んでもよい。

(付記３)
前記第１の光子数情報は、チェレンコフ光子数を更に含んでもよい。

（付記４）
前記特定部は、前記第１の発光確率モデルにより示される分布の累積分布に基づいて、前記第１のタイミングを特定してもよい。

（付記５）
前記特定部は、前記第１の発光確率モデルにより示される分布が極値を取る時刻に基づいて、前記第１のタイミングを特定してもよい。

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

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

（付記８）
前記算出部は、前記第１の光子数情報に基づいて、各時刻における発光の分布関数を生成し、生成した前記分布関数に基づいて、前記第１の発光確率モデルを算出してもよい。

（付記９）
前記取得部は、検出器毎に前記第１の光子数情報を取得してもよい。

（付記１０）
前記取得部は、イベントごとに前記第１の光子数情報を取得してもよい。

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

（付記１２）
本発明の一つの側面において提供されるデータ処理方法は、第１の検出器で検出される第１の光子数情報を取得し、前記第１の光子数情報に基づいて前記第１の検出器に対応する第１の発光確率モデルを算出し、前記第１の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第１のタイミングを特定し、前記第１の検出器で検出されるイベントの検出タイミングを測定し、前記第１のタイミングに基づいて、前記検出タイミングを補正する。

（付記１３）
本発明の一つの側面において提供されるプログラムは、第１の検出器で検出される第１の光子数情報を取得し、前記第１の光子数情報に基づいて前記第１の検出器に対応する第１の発光確率モデルを算出し、前記第１の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第１のタイミングを特定し、前記第１の検出器で検出されるイベントの検出タイミングを測定し、前記第１のタイミングに基づいて、前記検出タイミングを補正する処理をコンピュータに実行させる。

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

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

１５０ａ 取得機能
１５０ｂ 算出機能
１５０ｃ 特定機能
１５０ｄ 測定機能
１５０ｅ 補正機能
１５０ｆ 再構成機能
１５０ｇ 制御機能
１５０ｈ 受付機能
１５０ｉ 画像生成機能
１５０ｊ 表示制御機能

Claims (13)

1. 第１の検出器で検出される第１の光子数情報を取得する取得部と、
   前記第１の光子数情報に基づいて前記第１の検出器に対応する第１の発光確率モデルを算出する算出部と、
   前記第１の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第１のタイミングを特定する特定部と、
   前記第１の検出器で検出されるイベントの検出タイミングを測定する測定部と、
   前記第１のタイミングに基づいて、前記検出タイミングを補正する補正部と、
   を備える、核医学診断装置。
2. 前記第１の光子数情報は、シンチレーション光子数を含む、請求項１に記載の核医学診断装置。
3. 前記第１の光子数情報は、チェレンコフ光子数を更に含む、請求項２に記載の核医学診断装置。
4. 前記特定部は、前記第１の発光確率モデルにより示される分布の累積分布に基づいて、前記第１のタイミングを特定する、請求項１に記載の核医学診断装置。
5. 前記特定部は、前記第１の発光確率モデルにより示される分布が極値を取る時刻に基づいて、前記第１のタイミングを特定する、請求項１に記載の核医学診断装置。
6. 前記第１の発光確率モデルは、検出事象に含まれる複数の光子の内の最初の光子を検出する確率密度に関する、請求項１に記載の核医学診断装置。
7. 前記第１の発光確率モデルは、検出事象に含まれる複数の光子の内の所定の数の光子を検出する確率密度に関する、請求項１に記載の核医学診断装置。
8. 前記算出部は、前記第１の光子数情報に基づいて、各時刻における発光の分布関数を生成し、生成した前記分布関数に基づいて、前記第１の発光確率モデルを算出する、請求項６に記載の核医学診断装置。
9. 前記取得部は、検出器毎に前記第１の光子数情報を取得する、請求項１に記載の核医学診断装置。
10. 前記取得部は、イベントごとに前記第１の光子数情報を取得する、請求項１に記載の核医学診断装置。
11. 前記算出部は、シンチレーター内での発光位置の不定性を考慮して、前記第１の発光確率モデルを算出する、請求項１に記載の核医学診断装置。
12. 第１の検出器で検出される第１の光子数情報を取得し、
    前記第１の光子数情報に基づいて前記第１の検出器に対応する第１の発光確率モデルを算出し、
    前記第１の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第１のタイミングを特定し、
    前記第１の検出器で検出されるイベントの検出タイミングを測定し、
    前記第１のタイミングに基づいて、前記検出タイミングを補正する、データ処理方法。
13. 第１の検出器で検出される第１の光子数情報を取得し、
    前記第１の光子数情報に基づいて前記第１の検出器に対応する第１の発光確率モデルを算出し、
    前記第１の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第１のタイミングを特定し、
    前記第１の検出器で検出されるイベントの検出タイミングを測定し、
    前記第１のタイミングに基づいて、前記検出タイミングを補正する処理をコンピュータに実行させる、プログラム。

Images