(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024052403
(43)【公開日】2024-04-11
(54)【発明の名称】核医学診断装置、データ処理方法及びプログラム
(51)【国際特許分類】
G01T 1/161 20060101AFI20240404BHJP
G01T 1/22 20060101ALI20240404BHJP
【FI】
G01T1/161 E
G01T1/161 A
G01T1/22
【審査請求】未請求
【請求項の数】13
【出願形態】OL
(21)【出願番号】P 2022159094
(22)【出願日】2022-09-30
(71)【出願人】
【識別番号】594164542
【氏名又は名称】キヤノンメディカルシステムズ株式会社
(74)【代理人】
【識別番号】110001771
【氏名又は名称】弁理士法人虎ノ門知的財産事務所
(72)【発明者】
【氏名】河田 剛
(72)【発明者】
【氏名】勅使川原 学
【テーマコード(参考)】
2G188
4C188
【Fターム(参考)】
2G188AA02
2G188BB04
2G188CC23
4C188EE02
4C188FF04
4C188FF07
4C188GG17
4C188KK03
4C188KK13
4C188KK15
4C188KK35
4C188LL18
4C188LL30
(57)【要約】
【課題】画質を向上させること。
【解決手段】実施形態に係る核医学診断装置は、取得部と、算出部と、特定部と、測定部と、補正部とを備える。取得部は、第1の検出器で検出される第1の光子数情報を取得する。算出部は、前記第1の光子数情報に基づいて前記第1の検出器に対応する第1の発光確率モデルを算出する。特定部は、前記第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定する。測定部は、前記第1の検出器で検出されるイベントの検出タイミングを測定する。補正部は、前記第1のタイミングに基づいて、前記検出タイミングを補正する。
【選択図】
図1
【特許請求の範囲】
【請求項1】
第1の検出器で検出される第1の光子数情報を取得する取得部と、
前記第1の光子数情報に基づいて前記第1の検出器に対応する第1の発光確率モデルを算出する算出部と、
前記第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定する特定部と、
前記第1の検出器で検出されるイベントの検出タイミングを測定する測定部と、
前記第1のタイミングに基づいて、前記検出タイミングを補正する補正部と、
を備える、核医学診断装置。
【請求項2】
前記第1の光子数情報は、シンチレーション光子数を含む、請求項1に記載の核医学診断装置。
【請求項3】
前記第1の光子数情報は、チェレンコフ光子数を更に含む、請求項2に記載の核医学診断装置。
【請求項4】
前記特定部は、前記第1の発光確率モデルにより示される分布の累積分布に基づいて、前記第1のタイミングを特定する、請求項1に記載の核医学診断装置。
【請求項5】
前記特定部は、前記第1の発光確率モデルにより示される分布が極値を取る時刻に基づいて、前記第1のタイミングを特定する、請求項1に記載の核医学診断装置。
【請求項6】
前記第1の発光確率モデルは、検出事象に含まれる複数の光子の内の最初の光子を検出する確率密度に関する、請求項1に記載の核医学診断装置。
【請求項7】
前記第1の発光確率モデルは、検出事象に含まれる複数の光子の内の所定の数の光子を検出する確率密度に関する、請求項1に記載の核医学診断装置。
【請求項8】
前記算出部は、前記第1の光子数情報に基づいて、各時刻における発光の分布関数を生成し、生成した前記分布関数に基づいて、前記第1の発光確率モデルを算出する、請求項6に記載の核医学診断装置。
【請求項9】
前記取得部は、検出器毎に前記第1の光子数情報を取得する、請求項1に記載の核医学診断装置。
【請求項10】
前記取得部は、イベントごとに前記第1の光子数情報を取得する、請求項1に記載の核医学診断装置。
【請求項11】
前記算出部は、シンチレーター内での発光位置の不定性を考慮して、前記第1の発光確率モデルを算出する、請求項1に記載の核医学診断装置。
【請求項12】
第1の検出器で検出される第1の光子数情報を取得し、
前記第1の光子数情報に基づいて前記第1の検出器に対応する第1の発光確率モデルを算出し、
前記第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定し、
前記第1の検出器で検出されるイベントの検出タイミングを測定し、
前記第1のタイミングに基づいて、前記検出タイミングを補正する、データ処理方法。
【請求項13】
第1の検出器で検出される第1の光子数情報を取得し、
前記第1の光子数情報に基づいて前記第1の検出器に対応する第1の発光確率モデルを算出し、
前記第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定し、
前記第1の検出器で検出されるイベントの検出タイミングを測定し、
前記第1のタイミングに基づいて、前記検出タイミングを補正する処理をコンピュータに実行させる、プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本明細書及び図面に開示の実施形態は、核医学診断装置、データ処理方法及びプログラムに関する。
【背景技術】
【0002】
PET画像の再構成を行う際において、ToF(Time Of Flight)カーネルを用いて画像再構成が行われる。ToFカーネルとは、2つの検出器で検出された信号の検出時間差の関数として表された、検出イベントの確率密度関数である。PET画像の再構成には、精度が良いToFカーネルを用いて画像再構成を行うのが望ましい。また、PET画像の再構成において、ToFスペクトルの半値幅は、狭いほうが望ましい。
【0003】
このような観点から、ガンマ線がシンチレータに入射してから、発光が起こるまでの時間を精度よく推定できるのが望ましい。これに対して、例えば、検出される信号波形と検出時刻との間の関係をニューラルネットワークで学習する方法も考えられる。しかしながら、ニューラルネットワークを用いて学習を行う場合、検出される信号波形と検出時刻との間を機械学習する必要があるなど、多くの計算資源が必要となる場合もあった。
【先行技術文献】
【非特許文献】
【0004】
【非特許文献1】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年
【発明の概要】
【発明が解決しようとする課題】
【0005】
本明細書及び図面の開示の実施形態が解決しようとする課題の一つは、画質を向上させることである。ただし、本明細書及び図面に開示の実施形態により解決しようとする課題は上記課題に限られない。後述する実施形態に示す各構成による各効果に対応する課題を他の課題として位置づけることもできる。
【課題を解決するための手段】
【0006】
実施形態に係る核医学診断装置は、取得部と、算出部と、特定部と、測定部と、補正部とを備える。取得部は、第1の検出器で検出される第1の光子数情報を取得する。算出部は、前記第1の光子数情報に基づいて前記第1の検出器に対応する第1の発光確率モデルを算出する。特定部は、前記第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定する。測定部は、前記第1の検出器で検出されるイベントの検出タイミングを測定する。補正部は、前記第1のタイミングに基づいて、前記検出タイミングを補正する。
【図面の簡単な説明】
【0007】
【
図1】
図1は、実施形態に係る核医学診断装置の構成の一例を示した図である。
【
図2】
図2は、実施形態に係る核医学診断装置が行う処理の流れを示したフローチャートである。
【
図3】
図3は、
図2のステップS30の処理の一例をより詳細に説明したフローチャートである。
【
図4】
図4は、実施形態に係る核医学診断装置が行う処理について説明した図である。
【
図5】
図5は、
図3のステップS300の処理の一例をより詳細に説明したフローチャートである。
【
図6】
図6は、実施形態に係る核医学診断装置が行う処理について説明した図である。
【
図7】
図7は、実施形態に係る核医学診断装置が行う処理について説明した図である。
【
図8】
図8は、実施形態に係る核医学診断装置が行う処理について説明した図である。
【
図9】
図9は、
図2のステップS300の処理の別の一例をより詳細に説明したフローチャートである。
【
図10】
図10は、実施形態に係る核医学診断装置が行う処理について説明した図である。
【
図11】
図11は、実施形態に係る計算結果の一例について説明した図である。
【
図12】
図12は、実施形態に係る計算結果の一例について説明した図である。
【
図13】
図13は、実施形態に係る計算結果の一例について説明した図である。
【
図14】
図14は、実施形態に係る計算結果の一例について説明した図である。
【発明を実施するための形態】
【0008】
(実施形態)
以下、図面を参照しながら、核医学診断装置、データ処理方法及びプログラムの実施形態について詳細に説明する。
【0009】
図1は、実施形態に係る核医学診断装置としてのPET装置100の構成を示す図である。
図1に示すように、実施形態に係るPET装置100は、架台装置1と、コンソール装置2とを備える。架台装置1は、検出器3と、フロントエンド回路102と、天板103と、寝台104と、寝台駆動部106とを備える。
【0010】
検出器3は、被検体に投与され集積した薬剤から放出される陽電子が周囲の組織の電子と対消滅を起こすことで発生したガンマ線が発光体と相互作用することにより励起状態となった物質が再び基底状態に遷移する際に再放出される光であるシンチレーション光(蛍光)を検出することにより、放射線を検出する検出器である。また、実施形態においては、検出器3は、チェレンコフ光も検出することができる。検出器3は、被検体に投与され集積した薬剤から放出される陽電子が周囲の組織の電子と対消滅を起こすことで発生したガンマ線の放射線のエネルギー情報を検出する。検出器3は、被写体Pの周囲をリング状に取り囲むように複数配置され、例えば複数の検出器ブロックからなる。
【0011】
検出器3は、典型的には、シンチレータ結晶と、光検出素子からなる光検出面とからなる。
【0012】
シンチレータ結晶の材質としては、例えばチェレンコフ光を発生させるのに適した材質、例えば、ゲルマニウム酸ビスマス(BGO、Bismuth Germanium Oxide)や、鉛ガラス(SiO2+PbO)、フッ化鉛(PbF2)、PWO(PbWO4)等の鉛化合物を用いることができる。また、別の例として、例えば、LYSO(Lutetium Yttrium Oxyorthosilicate)、LSO(Lutetium Oxyorthosilicate)、LGSO(Lutetium Gadolinium Oxyorthosilicate)等やBGO等の、シンチレータ結晶を用いてもよい。光検出面3bを構成する光検出素子は、例えば複数の画素からなり、これらの画素それぞれは、例えばSPAD(Single Photon Avalanche Diode)で構成される。また、検出器3の構成は、上述の例に限られず、一例として、光検出素子は、例えばSiPM(Silicon photomultiplier)や、光電子増倍管を用いてもよい。
【0013】
なお、シンチレータ結晶はモノリシック結晶であってもよく、また光検出素子からなる光検出面は、シンチレータ結晶の6面に例えば配置されてもよい。なお、以下の実施形態では、検出器3のシンチレータ結晶は、モノリシック結晶でない場合をまずは例にとり説明する。
【0014】
また、架台装置1は、フロントエンド回路102により、検出器3の出力信号から計数情報を生成し、生成した計数情報を、コンソール装置2の記憶部130に格納する。なお、検出器3は、複数のブロックに区分けされ、フロントエンド回路102を備える。
【0015】
フロントエンド回路102は、検出器3の出力信号をデジタルデータに変換し、計数情報を生成する。この計数情報には、消滅ガンマ線の検出位置、エネルギー値、及び検出時間が含まれる。例えば、フロントエンド回路102は、シンチレーション光を同じタイミングで電気信号に変換した複数の光検出素子を特定する。そして、フロントエンド回路102は、消滅ガンマ線が入射したシンチレータの位置を示すシンチレータ番号(P)を特定する。消滅ガンマ線が入射したシンチレータ位置を特定する手段は、各光検出素子の位置及び電気信号の強度に基づいて重心演算を行うことによって特定してもよい。また、シンチレータと光検出素子の各々の素子サイズが対応している場合には、例えば最大出力が得られた光検出素子に対応するシンチレータを消滅ガンマ線が入射したシンチレータ位置と仮定し、さらにシンチレータ間散乱を考慮して最終的に特定するなどする。
【0016】
また、フロントエンド回路102は、各光検出素子から出力された電気信号の強度を積分計算するあるいは電気信号強度が閾値を超えた時間(Time over Threshold)を計測し検出器3に入射した消滅ガンマ線のエネルギー値(E)を特定する。また、フロントエンド回路102は、検出器3によって消滅ガンマ線によるシンチレーション光が検出された検出時間(T)を特定する。なお、検出時間(T)は、絶対時刻であってもよいし、撮影開始時点からの経過時間であってもよい。このように、フロントエンド回路102は、シンチレータ番号(P)、エネルギー値(E)、及び検出時間(T)を含む計数情報を生成する。
【0017】
なお、フロントエンド回路102は、例えば、CPU(Central Processing Unit)、GPU(Graphical Processing Unit)或いは、特定用途向け集積回路(Application Specific Integrated Circuit:ASIC)、プログラマブル論理デバイス(例えば、単純プログラマブル論理デバイス(Simple Programmable Logic Device:SPLD)、複合プログラマブル論理デバイス(Complex Programmable Logic Device:CPLD)、及びフィールドプログラマブルゲートアレイ(Field Programmable Gate Array:FPGA))等の回路により実現される。フロントエンド回路102は、フロントエンド部の一例である。
【0018】
天板103は、被写体Pが載置されるベッドであり、寝台104の上に配置される。寝台駆動部106は、処理回路150の制御機能105gによる制御の下、天板103を移動させる。例えば、寝台駆動部106は、天板103を移動させることで、被写体Pを架台装置1の撮影口内に移動させる。
【0019】
コンソール装置2は、操作者によるPET装置100の操作を受け付け、PET画像の撮影を制御するとともに、架台装置1によって収集された計数情報を用いてPET画像を再構成する。
図1に示すように、コンソール装置2は、処理回路150と、入力装置110と、ディスプレイ120と、記憶部130とを備える。なお、コンソール装置2が備える各部は、バスを介して接続される。処理回路150の詳細については後述する。
【0020】
入力装置110は、PET装置100の操作者によって各種指示や各種設定の入力に用いられるマウスやキーボード等であり、入力された各種指示や各種設定を、処理回路150に転送する。例えば、入力装置110は、撮影開始指示の入力に用いられる。
【0021】
ディスプレイ120は、操作者によって参照されるモニター等であり、処理回路150による制御の下、被写体の呼吸波形やPET画像を表示したり、操作者から各種指示や各種設定を受け付けるためのGUI(Graphical User Interface)を表示したりする。
【0022】
記憶部130は、PET装置100において用いられる各種データを記憶する。記憶部130は、例えば、メモリで構成され、一例として、RAM(Random Access Memory)、フラッシュメモリ(flash memory)等の半導体メモリ素子や、ハードディスク、光ディスク等によって実現される。記憶部130は、シンチレータ番号(P)、エネルギー値(E)、及び検出時間(T)が対応づけられた情報である計数情報、同時計数情報の通し番号であるコインシデンスNo.に計数情報の組が対応づけられた同時計数情報、または同時計数情報を集計して得られる射影データ、再構成されたPET画像等を記憶する。
【0023】
処理回路150は、取得機能150a、算出機能150b、特定機能150c、測定機能150d、補正機能150e、再構成機能150f、制御機能150g、受付機能150h、画像生成機能150i、表示制御機能150jを有する。なお、取得機能150a、算出機能150b、特定機能150c、補正機能150e及び再構成機能150fの各機能については、後ほど詳しく説明する。
【0024】
実施形態では、取得機能150a、算出機能150b、特定機能150c、測定機能150d、補正機能150e、再構成機能150f、制御機能150g、受付機能150h、画像生成機能150i、表示制御機能150jを有する。なお、取得機能150a、算出機能150b、特定機能150c、補正機能150e及び再構成機能150fにて行われる各処理機能は、コンピュータによって実行可能なプログラムの形態で記憶部130へ記憶されている。処理回路150はプログラムを記憶部130から読み出し、実行することで各プログラムに対応する機能を実現するプロセッサである。換言すると、各プログラムを読み出した状態の処理回路150は、
図1の処理回路150内に示された各機能を有することになる。
【0025】
なお、
図1においては単一の処理回路150にて、取得機能150a、算出機能150b、特定機能150c、測定機能150d、補正機能150e、再構成機能150f、制御機能150g、受付機能150h、画像生成機能150i、表示制御機能150jにて行われる処理機能が実現されるものとして説明するが、複数の独立したプロセッサを組み合わせて処理回路150を構成し、各プロセッサがプログラムを実行することにより機能を実現するものとしても構わない。換言すると、上述のそれぞれの機能がプログラムとして構成され、1つの処理回路150が各プログラムを実行する場合であってもよい。別の例として、特定の機能が専用の独立したプログラム実行回路に実装される場合であってもよい。
【0026】
上記説明において用いた「プロセッサ」という文言は、例えば、CPU(Central Processing Unit)、GPU(Graphical Processing Unit)或いは、特定用途向け集積回路(Application Specific Integrated Circuit:ASIC)、プログラマブル論理デバイス(例えば、単純プログラマブル論理デバイス(Simple Programmable Logic Device:SPLD)、複合プログラマブル論理デバイス(Complex Programmable Logic Device:CPLD)、及びフィールドプログラマブルゲートアレイ(Field Programmable Gate Array:FPGA))等の回路を意味する。プロセッサは記憶部130に保存されたプログラムを読み出し実行することで機能を実現する。
【0027】
なお、
図1において、取得機能150a、算出機能150b、特定機能150c、測定機能150d、補正機能150e、再構成機能150f、制御機能150g、受付機能150h、画像生成機能150i、表示制御機能150jは、それぞれ、取得部、算出部、測定部、補正部、再構成部、制御部、受付部、画像生成部、表示制御部の一例である。
【0028】
なお、処理回路150に代えて、フロントエンド回路102が、計測部、特定部等の処理を担ってもよい。
【0029】
処理回路150は、再構成機能150fにより、フロントエンド回路102から取得したデータに基づいて、PET画像の再構成を行い、画像生成機能150iにより、画像の生成を行う。
【0030】
処理回路150は、制御機能150gにより、架台装置1及びコンソール装置2を制御することによって、PET装置100の全体制御を行う。例えば、処理回路150は、制御機能150gにより、PET装置100における撮影を制御する。また、処理回路105は、制御機能150gにより、寝台駆動部106を制御する。
【0031】
処理回路150は、受付機能150hにより、入力装置110を通じて、ユーザから情報の入力を受け付ける。また、処理回路150は、表示制御機能150jにより、PET画像その他のデータをディスプレイ120に表示させる。また、処理回路150は、画像生成機能150iにより、種々の画像の生成を行う。
【0032】
続いて、実施形態に係る背景について簡単に説明する。
【0033】
PET画像の再構成を行う際において、ToF(Time Of Flight)カーネルを用いて画像再構成が行われる。ToFカーネルとは、2つの検出器で検出された信号の検出時間差の関数として表された、検出イベントの確率密度関数である。PET画像の再構成には、精度が良いToFカーネルを用いて画像再構成を行うのが望ましい。また、PET画像の再構成において、ToFスペクトルの半値幅は、狭いほうが望ましい。
【0034】
このような観点から、ガンマ線がシンチレータに入射してから、発光が起こるまでの時間を精度よく推定できるのが望ましい。これに対して、例えば、検出される信号波形と検出時刻との間の関係をニューラルネットワークで学習する方法も考えられる。しかしながら、ニューラルネットワークを用いて学習を行う場合、検出される信号波形と検出時刻との間を機械学習する必要があるなど、多くの計算資源が必要となる場合もあった。
【0035】
そこで、実施形態では、発光モデルに基づいて、ガンマ線がシンチレータに入射してから発光が起こるまでの遅延時刻を推定し、推定結果に基づいてToFカーネルの補正を行う。このことで、ToFスペクトルを先鋭化することが可能となり、再構成後のPET画像の画質を向上させることができる。
【0036】
具体的には、実施形態に係る核医学診断装置は、取得部と、算出部と、特定部と、測定部と、補正部とを備える。取得部は、第1の検出器で検出される第1の光子数情報を取得する。算出部は、第1の光子数情報に基づいて第1の検出器に対応する第1の発光確率モデルを算出する。特定部は、第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定する。測定部は、第1の検出器で検出されるイベントの検出タイミングを測定する。補正部は、第1のタイミングに基づいて、前記検出タイミングを補正する。
【0037】
また、実施形態に係るデータ処理方法は、第1の検出器で検出される第1の光子数情報を取得し、第1の光子数情報に基づいて第1の検出器に対応する第1の発光確率モデルを算出し、第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定し、第1の検出器で検出されるイベントの検出タイミングを測定し、第1のタイミングに基づいて、検出タイミングを補正する。
【0038】
また、実施形態に係るプログラムは、第1の検出器で検出される第1の光子数情報を取得し、第1の光子数情報に基づいて第1の検出器に対応する第1の発光確率モデルを算出し、第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定し、第1の検出器で検出されるイベントの検出タイミングを測定し、第1のタイミングに基づいて、検出タイミングを補正する処理をコンピュータに実行させる。
【0039】
以下、
図2~
図10を用いて、実施形態に係る核医学診断装置が行う処理について、詳しく説明する。
図2に、実施形態に係る核医学診断装置としてのPET装置100が行う処理の流れが示されている。以下のフローチャートでは、簡単のため、第1の検出器lと第2の検出器mに係る確率分布モデル(ToFカーネル)を生成する場合についてのみ記載している。しかしながら、典型的には、画像再構成においては、すべてのlとmの組について処理回路150は確率分布モデルを生成した上でステップS70の画像再構成処理を行う。この場合、処理回路150は、すべてのlとmの組についてステップS10からステップS60の処理を繰り返した上で、ステップS70の処理を行う。
【0040】
ステップS10において、処理回路150は、処理を行う検出器の組を示すlとmの値を選択する。すなわち、処理回路150は、第1のイベントを検出する検出器である第1の検出器l(l番目の検出器)と、第2のイベントを検出する検出器である第2の検出器m(m番目の検出器)とを特徴づけるlとmの値を選択する。なお、ここでいう第1のイベントと第2のイベントとで、一つの同時計数イベントが構成される。別の言い方をすれば、同時計数イベントの各々は、第1のイベントと第2のイベントの2つのイベントからなる。
【0041】
続いて、ステップS20において、処理回路150は、同時計数イベントの番号を示すiの値に1をセットする。ここで、i=1とは、処理回路150が、1番目の同時計数イベントの処理を行うことを意味する。
【0042】
続いて、ステップS30において、処理回路150は、算出機能150bにより、l番目の検出器で検出されたイベントと、m番目の検出器で検出されたイベントとの同時計数イベントのうちのi番目の同時計数イベントに対して、コインシデンスに関する確率分布モデルpct;l、m
i、(doi)(θl
i,θm
i、t)を算出する。ここで、θl
i、θm
iは、i番目の同時計数イベントにおいて、それぞれl番目の検出器及びm番目の検出器について取得されるパラメータを示す。また、(doi)の記号は、シンチレータ内での発光位置の不定性を考慮してもよいし、しなくてもよいことを意味し、シンチレータ内での発光位置の不定性を考慮する実施形態においては、処理回路150は、算出機能150bにより、シンチレータ内での発光位置の不定性を考慮したコインシデンスに関する確率分布モデルpct;l、m
i、doi(θl
i,θm
i、t)を算出し、シンチレータ内での発光位置の不定性を考慮しない実施形態においては、シンチレータ内での発光位置の不定性を考慮しないコインシデンスに関する確率分布モデルpct;l、m
i(θl
i,θm
i、t)を算出することを意味する。
【0043】
図3に、
図2のステップS30の処理について、より詳細な処理の流れが示されている。すなわち、
図3のステップS100~ステップS500は、
図2のステップS30をより詳細に説明したものとなる。
【0044】
はじめに、ステップS100において、処理回路150は、取得機能150aにより、i番目の同時計数イベントについて、第1の検出器lで検出される第1のイベントに関する第1の光子数情報を取得する。
【0045】
ここで、第1の光子数情報とは、例えば、当該イベントについてのシンチレーション光子数Ni
S;lである。また、別の例として、第1の光子数情報とは、例えば、当該イベントについてのチェレンコフ光子数Ni
C;lである。また、別の例として、第1の光子数情報とは、当該イベントについてのシンチレーション光子数Ni
S;l及びチェレンコフ光子数Ni
C;lである。すなわち、処理回路150は、取得機能150aにより、シンチレーション光子数Ni
S;l、チェレンコフ光子数Ni
C;l、またはシンチレーション光子数Ni
S;l及びチェレンコフ光子数Ni
C;lを、第1の光子数情報として取得する。
【0046】
なお、処理回路150が取得機能150aにより取得する第1の光子数情報(及び第2の光子数情報)は、典型的には、例えばイベントごとに取得されるものである。しかしながら、実施形態はこれに限られず、処理回路150は、取得機能150aにより、例えば第1の検出器lで検出される平均光子数等を、第1の光子数情報(及び第2の光子数情報)等として取得してもよい。
【0047】
また、処理回路150が取得機能150aにより取得する第1の光子数情報(及び第2の光子数情報)は、典型的には、例えば検出器ごとに取得されるものである。しかしながら、実施形態はこれに限られず、処理回路150は、取得機能150aにより、例えばすべての検出器にわたって同一の光子数情報を、第1の光子数情報(及び第2の光子数情報)等として取得してもよい。
【0048】
また、PET装置100は、第1の検出器lで取得されるエネルギースペクトルに基づいて、第1の光子数情報を取得してもよい。
【0049】
図4に、検出器3内での発光モデルと、処理回路150が取得機能150aにより取得する光子スペクトルに関する説明が示されている。検出器3のシンチレータ11にガンマ線10が入射すると、ガンマ線10はシンチレータ11にエネルギーを光電吸収またはコンプトン散乱の形で付与するとともに、高エネルギー電子12が発生する。曲線16は、高エネルギー電子12のエネルギーE
eを時間tの関数として示したものである。高エネルギー電子12は、シンチレータ11中を走行するとともにエネルギーを失っていくが、エネルギーが閾値E
thより大きい間は、チェレンコフ光13を周囲に放出する。
【0050】
曲線17は、チェレンコフ光の光子密度nc(t)を、時間tのガウス関数に近似したものである。具体的には、チェレンコフ光の光子密度nc(t)は、Ncをチェレンコフ光の光子数、σCをチェレンコフ光スペクトルの半値幅、t0をチェレンコフ光のピーク時刻として、以下の式(1)で与えられる。
【0051】
【0052】
一方、高エネルギー電子12は、走行に伴い、シンチレータ11の中で価電子14を励起し、励起された価電子14が基底状態に戻るとき、シンチレーション光15が放出される。曲線19は、価電子の励起密度nex(t)を、時刻tのガウス関数で近似したしたものである。価電子の励起密度nex(t)は、Nsをシンチレーション光の光子数、σSをシンチレーション光スペクトルの半値幅、t1をシンチレーション光のピーク時刻として、以下の式(2)で与えられる。
【0053】
【0054】
曲線21は、時刻20において励起された価電子がもとの状態に遷移していく様子を示している。励起された価電子がもとの状態に遷移するのに伴ってシンチレーション光が発生する。曲線22は、シンチレーション光の光子密度ns(t)を、時刻tの関数として示したものである。シンチレーション光の光子密度ns(t)は、緩和定数をτSとして、価電子の励起密度nex(t)を用いて、以下の式(3)で与えられる。
【0055】
【0056】
ここで、価電子の励起密度nex(t)の関数形を式(2)で仮定した場合、式(3)に式(2)を代入すると、以下の式(4)が得られる。
【0057】
【0058】
すなわち、一例として、PET装置100は、第1の検出器l及び第2の検出器mの各々で取得されるエネルギースペクトルを取得し、当該エネルギースペクトルを基に曲線17や曲線22に示されるスペクトルを生成し、生成したスペクトルに基づいて、第1の光子数情報及び第2の光子数情報を取得してもよい。
【0059】
なお、発光の分布関数pph(t)は、チェレンコフ光の光子密度nc(t)とシンチレーション光の光子密度ns(t)との和に比例するから、以下の式(5)が成り立つ。
【0060】
【0061】
ここで、発光の分布関数pph(t)は、以下の式(6)のように規格化されている。
【0062】
【0063】
なお、ステップS100において、処理回路150は、取得機能150aにより、光子数情報以外の情報を、第1の検出器lに係るパラメータセットとして、取得してもよい。一例として、処理回路150は、取得機能150aにより、チェレンコフ光に対する光センサの時間分解能σCを、第1の光子数情報に代えて、または第1の光子数情報に加えて、第1の検出器lにおいて取得するパラメータセットθlの一つとして、取得してもよい。また、別の例として、処理回路150は、取得機能150aにより、シンチレーション光に対する光センサの時間分解能σSを、第1の光子数情報に代えて、または第1の光子数情報に加えて、第1の検出器lにおいて取得するパラメータセットθlの一つとして、取得してもよい。また、別の例として、処理回路150は、取得機能150aにより、過剰ノイズ発生確率を、第1の検出器lにおいて取得するパラメータセットθlの一つとして、取得してもよい。ここでいう過剰ノイズの発生確率とは、例えばオプティカルクロストークやアフターパルス等のAPD(Avalanche photo diode)過剰ノイズの発生確率を言う。
【0064】
ステップS100と同様に、ステップS200において、処理回路150は、取得機能150aにより、i番目の同時計数イベントについて、第1の検出器lとは異なる第2の検出器mで検出される第2のイベントに関する第2の光子数情報を取得する。
【0065】
ここで、第2の光子数情報とは、例えば、当該同時計数イベントについてのシンチレーション光子数Ni
S;mである。また、別の例として、第2の光子数情報とは、例えば、当該同時計数イベントについてのチェレンコフ光子数Ni
C;mである。また、別の例として、第2の光子数情報とは、当該同時計数イベントについてのシンチレーション光子数Ni
S;m及びチェレンコフ光子数Ni
C;mである。すなわち、処理回路150は、取得機能150aにより、シンチレーション光子数Ni
S;m、チェレンコフ光子数Ni
C;m、またはシンチレーション光子数Ni
S;m及びチェレンコフ光子数Ni
C;mを、第2の光子数情報として取得する。ここで、処理回路150が取得機能150aにより計測する第2の光子数情報は、典型的には、イベントごとに取得されるものである。
【0066】
PET装置100は、第2の検出器mで取得されるエネルギースペクトルに基づいて、第2の光子数情報を取得してもよい。
【0067】
また、処理回路150は、取得機能150aにより、チェレンコフ光に対する光センサの時間分解能σC、シンチレーション光に対する光センサの時間分解能σS、過剰ノイズ発生確率等を、第2の検出器mにおいて取得するパラメータセットθmとして、取得してもよい。
【0068】
ステップS100に続き、ステップS300において、処理回路150は、算出機能150bにより、ステップS100で取得された第1の光子数情報に基づいて、i番目の同時計数イベントにおける、第1の検出器lに対する第1の発光確率モデルpt;l
i,(doi)を算出する。なお、処理回路150は、算出機能150bにより、第1の検出器lに対する第1の発光確率モデルpt;l
i,(doi)を、典型的には、同時計数イベントごとに定める。以下の実施形態では、処理回路150は、算出機能150bにより、第1の検出器lに対する第1の発光確率モデルpt;l
i,(doi)を、同時計数イベントごとに定めるとして説明する。しかしながら、実施形態はこれに限られず、処理回路150は、算出機能150bにより、第1の検出器lに対する第1の発光確率モデルを検出器単位でのみ定めてもよい。
【0069】
図5にステップS300の処理の一例がより詳細に説明されている。
図5のフローチャートは、ステップS300の処理の一例をより詳細に説明したフローチャートである。
【0070】
ステップS310において、処理回路150は、算出機能150bにより、ステップS100で取得された第1の検出器lにおける第1の光子数情報に基づいて、各時刻における発光の分布関数モデルpi
ph,l (t)を生成する。発光の分布関数pph (t)は、時刻tに発光が起こる確率を示す分布関数である。発光の分布関数モデルpi
ph,l (t)は、時刻に関して全時間で積分すると1になる。また、発光の分布関数モデルpi
ph,l(t)は、検出器ごと及び同時計数イベントごとに定められる。
【0071】
処理回路150は、算出機能150bにより、第1の光子数情報に基づいて、例えば式(1)、(4)、(5)、(6)等を用いて発光の分布関数モデルp
i
ph,l (t)を生成する。第1の光子数情報が、検出器ごと及び同時計数イベントごとに取得されているので、処理回路150は、算出機能150bにより、発光の分布関数モデルp
ph,l
i (t)を、検出器ごと及び同時計数イベントごとに算出する。
図6に、発光の分布関数p
ph(t)の一例が示されている。
図6において、曲線30は、発光の分布関数p
ph(t)を示している。また、検出時刻を定義するためのフォトンの閾値(N’)を併記する。
【0072】
続いて、ステップS320Aにおいて、処理回路150は、算出機能150bにより、発光の分布関数モデルp
i
ph,l (θ
i
l,t)に基づいて、第1の発光確率モデルp
i
t,l(≧N’,θ
i
l, t)を算出する。ここで、第1の発光確率モデルp
i
t,l(≧N’,θ
i
l, t)は、検出事象に含まれる複数の光子の内の最初のN’個以上の光子を検出する確率密度に関するモデルである。換言すると、第1の発光確率モデルp
t,l(≧N’,t)は、検出事象に含まれる複数の光子の内の所定の数の光子を検出する確率密度に関するモデルである。
図7に、第1の発光確率モデルp
t,l(Δt)が示されている。ここで、光子の検出閾値N’はある値に固定しているため、図示していない。
図7の曲線31は、第1の発光確率モデルp
t,l(Δt)の概形を示している。なお、第1の発光確率モデルp
t,l(Δt)は、検出器ごと及び同時計数イベントごとに定められる。以下、必要に応じて、検出器の添え字、同時計数イベントの番号を表す添え字、及び取得されるパラメータの集合を表す引数θ
i
l等については適宜省略する。なお、曲線32は、第1の発光確率モデルp
t,l(t)の累積分布P(t)を表す曲線である。曲線32については後述する。
【0073】
続いて、第1の発光確率モデルpt,l(≧N’,t)の算出方法について説明する。検出事象に含まれる複数の光子の内の最初のN’個以上の光子を検出する確率密度である第1の発光確率モデルpt(≧N’,t)は、以下のように算出することができる。
【0074】
はじめに、時刻tまでにN’個未満のフォトンしか検出しない確率P’(≧N’,t)に関する方程式を導出する。これは、時刻tまでに高々N’-1個のフォトンしか検出していない確率であり、時刻tまでにk個のフォトンを検出している確率Pd(k,t)を用いて以下の式(7)で表される。
【0075】
【0076】
時刻tまでにk個のフォトンを検出している確率は、フォトンの検出事象を2項分布で評価することで、以下の式(8)のように微分形式で表現できる。
【0077】
【0078】
この微分方程式を順次k=0、1,2,…に対して解くと、帰納的に以下の式(9)が成り立つ。
【0079】
【0080】
式(7)より、時刻tまでにN’個未満のフォトンしか検出しない確率P’(≧N’,t)は以下の式(10)のように書くことが出来る。
【0081】
【0082】
これを式(8)に代入することで、検出事象に含まれる複数の光子の内の最初のN’個の光子を検出する確率密度である第1の発光確率モデルpt(≧N’,t)は、以下の式(11)のように表すことができる。
【0083】
【0084】
ここで、関数F(t)は発光の分布関数pph(t)の累積分布関数であり、以下の式(12)のように書くことが出来る。
【0085】
【0086】
以下の実施例では簡単のため、初めてフォトンを検出した時刻で評価した確率分布、すなわち、N’=0とした場合を例にとる。すなわち、初の光子を検出する確率密度である第1の発光確率モデルをpt (t)とし、以下の式(13)のように陽に表す。
【0087】
【0088】
すなわち、処理回路150は、算出機能150bにより、式(12)及び式(13)を用いることで、発光の分布関数pph(t)に基づいて、第1の発光確率モデルpt(t)を計算することができる。このように、処理回路150は、算出機能150bにより、式(12)及び式(13)に第1の検出器における発光の分布関数モデルpph (t)を代入し、第1の光子数情報を式(13)の右辺のNに代入して式(13)の右辺を評価することで、第1の光子数情報及び同時計数イベントごとに定められた第1の検出器における発光の分布関数モデルpph (t)に基づいて、第1のイベントについての第1の発光確率モデルpt(t)を算出する。
【0089】
なお、上述の発光の分布関数pph(t)及び第1の第1の発光確率モデルpt(t)は、各検出器ごと及び同時計数イベントごとに定めることができる。従って、検出器の添え字及び同時計数イベントに関する添え字を明示的に記載すると、式(13)は、以下の式(14)のように書くこともできる。
【0090】
【0091】
ここで、pi
t;l(t)はi番目の同時計数イベントの検出器lにおける第1の発光確率モデル、pi
ph;l(t)はi番目の同時計数イベントの検出器lにおける発光の分布関数、Ni
s;lはi番目の同時計数イベントの検出器lにおけるシンチレーション光の光子数、Ni
c;lはi番目の同時計数イベントの検出器lにおけるチェレンコフ光の光子数を示す。
【0092】
なお、第1の発光確率モデルを算出の際、処理回路150は、算出機能150bにより、シンチレータ内での発光位置の不定性を考慮した計算を行うこともできる。例えば、
図8に示されるように、ガンマ線40がシンチレータ11の深さxにある点42で発光し、矢印43で示されるように、シンチレーション光が物質中の光速c/nで屈折率nの媒体中を伝播してシンチレータ11の端44に到達した場合、シンチレータ内の発光位置によって、ガンマ線40がシンチレータに入射してからシンチレータ11の端に到達するまでの時刻が異なるので、信号を検出する時刻が異なる。
【0093】
すなわち、ガンマ線40がシンチレータ11に入射してから深さxの地点まで到達するのに要する時刻t1は、以下のように、式(15)で与えられる。
【0094】
【0095】
また、点42で発光した光が深さxからシンチレータ11の端まで到達するのに要する時刻t2は、以下のように、式(16)で与えられる。
【0096】
【0097】
従って、ガンマ線40がシンチレータ11に入射してからシンチレーション光がシンチレータ11の端まで到達するのに要する時刻t’は、以下のように、式(17)で与えられる。
【0098】
【0099】
式(17)より、シンチレータ内の発光位置xによって、信号を検出する時刻が異なることがわかる。
【0100】
以上をもとに、式(13)を基に、シンチレータ内の発光位置xを考慮した定式化を考えると、シンチレータ内の発光位置xを考慮した取り扱いを考える場合、以下の式(18)が成り立つ。
【0101】
【0102】
ここで、pt、l
i,doi(t)は、同時計数イベントごとに定められた、シンチレータ内での発光位置の不定性を考慮した第1の発光確率モデルであり、pph、l
i,doi(t)は、同時計数イベントごとに定められた、シンチレータの位置の発光位置の不定性を考慮した発光の分布関数モデルである。
【0103】
なお、式(18)は、処理回路150が取得機能150aによりステップS100において取得するパラメータθlを関数の引数として明示的に表現すると、以下の式(19)及び式(20)のようにも書くことができる。
【0104】
【0105】
【0106】
ここで、θi
lは、同時計数イベントごとに定められた、検出器lにおいて取得されるパラメータセットを示しており、例えば、同時計数イベントごとに定められたシンチレーション光子数Ni
s;l、チェレンコフ光子数Ni
c;l、チェレンコフ光に対する光センサの時間分解能σi
C、シンチレーション光に対する光センサの時間分解能σi
S等のパラメータの組を表す。
【0107】
式(18)に戻り、シンチレータの位置の発光位置の不定性を考慮した場合のステップS300の処理について説明する。
図9は、シンチレータ内での発光位置の不定性を考慮した場合の、発光確率モデルの算出の流れを示したフローチャートである。
【0108】
図5で示した場合と同様に、ステップS310において、処理回路150は、算出機能150bにより、第1の検出器lにおける第1の光子数情報等のパラメータθ
lに基づいて、発光の分布関数モデルp
ph,l(t)を生成する。
【0109】
続いて、ステップS311において、処理回路150は、算出機能150bにより、シンチレータ内での位置xごとに、発光確率密度pdoi(x)を、減弱係数μの値に基づいて算出する。位置xにおける発光確率密度pdoi(x)は、以下の式(21)で表される。
【0110】
【0111】
なお、発光確率密度pdoi(x)は、以下の式(22)に示されるように、規格化されている。
【0112】
【0113】
なお、
図8において、グラフ45には、減弱率μexp(-μx)が示されている。
【0114】
図9に戻り、ステップS312において、処理回路150は、算出機能150bにより、ステップS310で生成された、同時計数イベントごとに定められた発光の分布関数である分布関数モデルp
i
ph、l(t)と発光確率密度p
doi(x)とに基づいて、同時計数イベントごとに定められた発光位置不確定性を考慮した発光の分布関数モデルp
ph;l
i,doi(t)を第1の発光確率モデルとして算出する。
【0115】
具体的には、処理回路150は、算出機能150bにより、発光位置不確定性を考慮した発光の分布関数モデルpi,doi
ph;l(t)を、ステップS310で生成された分布関数モデルpi
ph、l(t)と、ステップS311で算出された発光確率密度pdoi(x)とに基づいて、次の式(23)に基づいて算出する。
【0116】
【0117】
ここで、Lはシンチレータの長さであり、xは発光位置であり、t’は式(17)において説明した、ガンマ線がシンチレータに入射してからシンチレーション光がシンチレータの端まで到達するまでの遅延時間である。
【0118】
続いて、ステップS320Bにおいて、処理回路150は、算出機能150bにより、発光位置不確定性を考慮した発光の分布関数モデルpi,doi
ph;l(t)に基づいて、同発光位置不確定性を考慮した発光確率モデルpi,doi
t;l(θl;t)を算出する。具体的には、式(24)及び式(25)に示されるように、処理回路150は、算出機能150bにより、式(13)と同様に、発光位置不確定性を考慮した発光の分布関数モデルpi,doi
ph(t)に基づいて、発光位置不確定性を考慮した発光確率モデルpi,doi
t(t)を算出する。
【0119】
【0120】
【0121】
なお、式(23)において、式(18)を基に被積分関数のxをt’に変数変換すると、以下の式(26)が得られるので、処理回路150は、算出機能150bにより、式(25)の右辺の評価において、式(26)を用いて算出を行ってもよい。
【0122】
【0123】
このようにして、処理回路150は、算出機能150bにより、シンチレータ内での発光位置の不定性を考慮して、第1の発光確率モデルおよび第2の発光確率モデルを算出する。
【0124】
図3に戻り、ステップS400において、処理回路150は、算出機能150bにより、ステップS300と同様の処理を行うことにより、ステップS200で取得された第2の光子数情報に基づいて、第2の検出器mに対する、同時計数イベントごとに定められた第2の発光確率モデルp
i、(doi)
t;mを算出する。
【0125】
具体的には、処理回路150は、算出機能150bにより、ステップS200で取得された第2の検出器mにおける同時計数イベントごとに定められた第2の光子数情報に基づいて、例えば式(1)、(4)、(5)、(6)等に相当する式を用いて、同時計数イベントごとに定められた発光の分布関数モデルpi
ph;m (t)を生成する。
【0126】
続いて、処理回路150は、算出機能150bにより、第2の光子数情報及び第2の検出器mにおける発光の分布関数モデルpi
ph;m(t)に基づいて、同時計数イベントごとに定められた第2の発光確率モデルpi
t;m(t)を、例えば式(14)等に相当する式を用いて算出する。
【0127】
また、第2の発光確率モデルを算出の際、処理回路150は、算出機能150bにより、シンチレータ内での発光位置の不定性を考慮した計算を、式(17)及び式(21)~(26)と同様の式を用いて行うこともできる。この場合、処理回路150は、算出機能150bにより、シンチレータ内での位置xごとに、発光確率密度pdoi(x)を、減弱係数μの値に基づいて算出する。続いて、処理回路150は、算出機能150bにより、同時計数イベントごとに定められた分布関数モデルpi
ph;m(t)と発光確率密度pdoi(x)とに基づいて、同時計数イベントごとに定められた発光位置不確定性を考慮した発光の分布関数モデルpi、doi
ph;m(t)を算出する。続いて、処理回路150は、算出機能150bにより、同時計数イベントごとに定められた発光位置不確定性を考慮した発光の分布関数モデルpi、doi
ph;m(t)に基づいて、同時計数イベントごとに定められた発光位置不確定性を考慮した発光確率モデルpi、doi
t;m(t)を算出する。
【0128】
このようにして、処理回路150は、算出機能150bにより、シンチレータ内での発光位置の不定性を考慮して、同時計数イベントごとに定められた第2の発光確率モデルpt,m(t)を算出する。
【0129】
続いて、ステップS410において、処理回路150は、特定機能150cにより、第1の発光確率モデルに基づいて、発光の検出確率が所定の閾値以上となる第1のタイミングを特定する。一例として、処理回路150は、特定機能150cにより、ステップS300において算出された第1の発光確率モデルに基づいて、処理回路150は、特定機能150cにより、イベントごとの推定検出遅延時間Δtl、i
pを、第1のタイミングとして特定する。
【0130】
ここで、推定検出遅延時間とは、例えばX線がシンチレータに入射してから、発光の信号が検出されるまでにかかる時刻の推定値である。ここで、処理回路150は、特定機能150cにより、検出器ごと、同時計数イベントごとに、推定検出遅延時間の推定を行うことができる。処理回路150は、推定検出遅延時間を推定して検出時刻から差し引くことで、ToFスペクトルの先鋭化を実現することができ、画質を向上させることができる。
【0131】
第1のタイミングの特定について、
図7に戻り説明する。発光の検出確率が所定の閾値以上となる第1のタイミングを特定する第1の方法として、処理回路150は、特定機能150cにより、第1の発光確率モデルp
i,doi
t;l(θ
l;t)により示される分布の累積分布P(t)に基づいて、発光の検出確率が所定の閾値以上となる第1のタイミングを特定する。ここで、第1の発光確率モデルp
i,doi
t;l(θ
l;t)により示される分布の累積分布P(t)は、例えば、以下の式(27)で与えられる。
【0132】
【0133】
ここで、処理回路150は、特定機能150cにより、第1の発光確率モデルpi,doi
t;l(θl;t)により示される分布の累積分布P(t)が、所定の閾値pに達するような時刻、すなわち曲線32と直線33との交点となる時刻34を、第1のタイミングΔtp
l,iとして特定する。これを式で表現すると、第1のタイミングΔtp
l,iは、以下の式(28)で与えられる。
【0134】
【0135】
また、発光の検出確率が所定の閾値以上となる第1のタイミングを特定する第2の方法として、処理回路150は、特定機能150cにより、第1の発光確率モデルpi,doi
t;l(θl;t)により示される分布が極値を取る時刻35を、第1のタイミングΔtp
l,iとして特定する。これを式で表現すると、第1のタイミングΔtp
l,iは、以下の式(29)で与えられる。
【0136】
【0137】
処理回路150は、特定機能150cにより、第1のタイミングΔtp
l,iを、典型的には、同時計数イベントごと及び検出器ごとに特定する。しかしながら、実施形態はこれに限られず、処理回路150は、特定機能150cにより、イベント単位以外で、または検出器単位以外で、発光の検出確率が所定の閾値以上となる第1のタイミングを特定してもよい。
【0138】
図3に戻り、ステップS410と同様にして、ステップS420において、処理回路150は、特定機能150cにより、第2の発光確率モデルに基づいて、発光の検出確率が所定の閾値以上となる第2のタイミングを特定する。一例として、処理回路150は、特定機能150cにより、ステップS300において算出された第2の発光確率モデルに基づいて、処理回路150は、特定機能150cにより、イベントごとの推定検出遅延時間Δt
m、i
pを、第2のタイミングとして特定する。
【0139】
発光の検出確率が所定の閾値以上となる第2のタイミングを特定する第1の方法として、処理回路150は、特定機能150cにより、第2の発光確率モデルpi,doi
t;m(θm;t)により示される分布の累積分布P(t)に基づいて、発光の検出確率が所定の閾値以上となる第1のタイミングを特定する。一例として、処理回路150は、特定機能150cにより、第1の発光確率モデルpi,doi
t;m(θm;t)により示される分布の累積分布P(t)が、所定の閾値pに達するような時刻を、第2のタイミングΔtp
m,iとして特定する。
【0140】
また、発光の検出確率が所定の閾値以上となる第2のタイミングを特定する第2の方法として、処理回路150は、特定機能150cにより、第2の発光確率モデルpi,doi
t;m(θm;t)により示される分布が極値を取る時刻を、第2のタイミングΔtp
m,iとして特定する。
【0141】
ステップS450において、処理回路150は、補正機能150eにより、ステップS410で推定された第1のタイミングΔtp
l,iに基づいて、検出器lで検出されたイベントであってi番目の同時計数イベントと対応付けられているイベントについて、検出時刻の補正を行う。ここで、補正後の検出時刻tl,i
corrected は、補正前の検出時刻tl,i
dと、ステップS410で特定された第1のタイミングΔtp
l,iとを用いて、以下の式(30)で表される。
【0142】
【0143】
換言すると、処理回路150は、測定機能150dにより、第1の検出器lで検出されるイベントの検出タイミングtl,i
dを測定し、補正機能150eに基づいて、第1のタイミングΔtp
l,iに基づいて、検出タイミングtl,i
dを補正し、補正後の検出時刻tl,i
correctedを得る。
【0144】
また、これに伴い、処理回路150は、補正機能150eにより、検出時刻が補正された第1の発光確率モデルpi,doi
t;l;corrected(θm;t)を算出する。ここで、補正された第1の発光確率モデルpi,doi
t;l;corrected(θl;t)と、補正前の第1の発光確率モデルpi,doi
t;l;corrected(θl;t)との関係は、例えば以下の式(31)で与えられる。
【0145】
【0146】
なお、説明を簡単にするために、ステップS500以降の処理において、補正された第1の発光確率モデルpi,doi
t;l;corrected(θl;t)を、単に、第1の発光確率モデルpi,doi
t;l(θl;t)と記載するが、以降のコインシデンスに関する確率分布モデルを算出する処理において、第1の発光確率モデルpi,doi
t;l(θl;t)は、補正された第1の発光確率モデルpi,doi
t;l;corrected(θl;t)を意味する。
【0147】
ステップS460において、ステップS450と同様の処理を行い、処理回路150は、補正機能150eにより、ステップS400で推定された第2のタイミングΔtp
m,iに基づいて、検出器mで検出されたイベントであってi番目の同時計数イベントと対応付けられているイベントについて、検出時刻の補正を行う。すなわち、処理回路150は、測定機能150dにより、第1の検出器lで検出されるイベントの検出タイミングtm,i
dを測定し、補正機能150eに基づいて、第2のタイミングΔtp
m,iに基づいて、検出タイミングtm,i
dを補正し、補正後の検出時刻tm,i
correctedを得る。
【0148】
また、これに伴い、処理回路150は、補正機能150eにより、検出時刻が補正された第2の発光確率モデルpi,doi
t;m;corrected(θm;t)を算出する。ステップS500以降の処理において、補正された第2の発光確率モデルpi,doi
t;m;corrected(θm;t)を、単に、第1の発光確率モデルpi,doi
t;m(θm;t)と記載するが、ステップS500以降のコインシデンスに関する確率分布モデルを算出する処理において、第2の発光確率モデルpi,doi
t;m(θm;t)は、補正された第2の発光確率モデルpi,doi
t;m;corrected(θm;t)を意味する。
【0149】
続いて、ステップS500において、処理回路150は、算出機能150bにより、ステップS300で算出された第1の発光確率モデルpt;l
doi(t1)と算出された第2の発光確率モデルpt;m
doi(t2)とに基づいて、同時計数イベントごとに特定された、第1のイベントと第2のイベントとに基づいて規定されるコインシデンスに関する確率分布モデルpi、(doi)
ct;l,m(θl, θm,t)を特定する。
【0150】
図10において、曲線53は、第1の検出器lのシンチレータ51における同時計数イベントごとに定められた第1の発光確率モデルp
i,doi
t;l(θ
l,t
l)を、遅延時間t
lの関数としてプロットしたものを示している。また、曲線54は、第2の検出器mのシンチレータ52における第2の発光確率モデルp
i,doi
t;m(θ
m,t
m)を、遅延時間t
mの関数としてプロットしたものを示している。処理回路150は、算出機能150bにより、検出器lにおける検出時刻t
lにおける事象と、検出器mにおける検出時刻t
mにおける積事象を、検出時間差の関数として表現した、同時計数イベントごとに定められた確率分布モデルp
i,doi
ct;l,m(θ
l, θ
m,t)を、特定する。曲線56は、このような確率分布モデルp
i,doi
ct;l,m(θ
l, θ
m,t)の例を示している。
【0151】
ここで、時間幅dtの間における、検出器lにおける検出時刻tlにおける事象と、検出器mにおける検出時刻tmにおける積事象の数は、以下の式(32)で与えられる。
【0152】
【0153】
ここで、θi
lは、i番目の同時計数イベントに係る、検出器lにおいて取得されるパラメータセットを示している。また、パラメータθを、より具体的な形で表現すると、以下の式(33)が成り立つ。
【0154】
【0155】
ここで、Ni
s;l、Ni
s;mは、i番目の同時計数イベントに係る、検出器l及びmにおけるシンチレーション光子数であり、Ni
c;l、Ni
c;mは、i番目の同時計数イベントに係る、検出器l及びmにおけるチェレンコフ光子数を示す。
【0156】
また、検出器lにおける検出時刻tlと、検出器mにおける検出時刻tmとの検出時間の差は、以下の式(34)で与えられる。
【0157】
【0158】
従って、pi,doi
ct;l,m(θl, θm,t)は、以下の式(35)で与えられる。
【0159】
【0160】
また、パラメータθを、より具体的な形で表現すると、以下の式(36)が成り立つ。
【0161】
【0162】
図2に戻り、ステップS30の処理により、処理回路150は、算出機能150bにより、第1の検出器lと第2の検出器mとによって規定されるLORの検出時間差に関する確率分布モデルp
i、doi
ct;l,m(θ
i
l, θ
i
m,t)を、i番目の同時計数イベントについて算出する。処理回路150は、この処理を、全ての同時計数イベントについて算出し、それらを加算することで、最終的な、第1の検出器lと第2の検出器mとによって規定されるLORの検出時間差に関する確率分布モデルp
doi
ct;l,m(θ
i
l, θ
i
m,t)を算出し、それを基に、データの再構成を行うことができる。
【0163】
ステップS40において、処理回路150は、算出機能150bにより、すべての同時計数イベントについて、ステップS30の処理が終わったか否かを判定する。処理回路150が、算出機能150bにより、すべての同時計数イベントについて、ステップS30の処理が終わっていないと判定した場合(ステップS40 No)、処理はステップS50に進み、i+1番目のデータについて、処理が行われる。一方、処理回路150が、算出機能150bにより、すべての同時計数イベントについて、ステップS30の処理が終わったと判定した場合(ステップS40 Yes)、処理はステップS60に進む。
【0164】
ステップS60において、処理回路150は、算出機能150bにより、ステップS30で算出された、各同時計数イベントごとの確率分布モデルpi、doi
ct;l,m(θi
l, θi
m,t)を加算して、第1の検出器lと第2の検出器mにおけるコインシデンスに関する確率分布モデルptot
ct;l、mを生成する。
【0165】
具体的には、第1の検出器lと第2の検出器mにおけるコインシデンスに関する確率分布モデルptot
ct;l、mは、以下の式(37)で表される。
【0166】
【0167】
続いて、ステップS70において、処理回路150は、再構成機能150fにより、ステップS60で特定された、第1の光子数情報と第2の光子数情報とに基づいて規定されたコインシデンスに関する確率分布モデルptot
ct;l、mをTOFカーネルとして利用することで、第1の検出器lと第2の検出器mとの間のLORに対応するデータの再構成を行う。
【0168】
図11~
図14に、実施形態に係る手法の検証結果が示されている。この検証では、第1のタイミングΔt
p
l,iを算出するための閾値pの値をさまざまに変化させ、再構成画像において放射線物質の分布のピークが先鋭化するかどうかを確認した。
【0169】
図11に、事象1、2、3の3つの事象について、第1の発光確率モデルp
i,doi
t;l(θ
l;t)が示されている。この3つの事象について、式(28)の右辺の閾値pを変化させて、それぞれの事象についての第1のタイミングΔt
pを算出した。直線60がp=0.1に対応し、直線61がp=0.4に対応し、直線62がp=0.8に対応する。また、矢印63は、p=0.4の時の第1の事象に対する第1のタイミングΔt
p_1を、矢印64は、p=0.4の時の第2の事象に対する第1のタイミングΔt
p_2を、矢印65は、p=0.4の時の第3の事象に対する第1のタイミングΔt
p_3を示している。
【0170】
図12に、p=0.4及びp=0.1の場合における、補正後のコインシデンススペクトルが示されている。グラフ71はp=0.4の時の補正後のコインシデンススペクトルであり、グラフ70はp=0.1の時の補正後のコインシデンススペクトルを示す。また、遅延時間が0付近でのカウント数を比較すると、p=0.1においてはカウント数が多いが幅も大きく、これに対して、p=0.4においては、カウント数も多く幅も小さくなった。従って、閾値をp=0.4付近にした場合、コインシデンススペクトルにおいてカウント数が大きくなりまた幅も小さくなる。
【0171】
また、実施形態に係る方法による補正を行った場合と行わない場合とで、再構成画像において、放射線物質の分布のピークが先鋭化するか否かを検証した。まず、
図13に示されるように、実測データ80に対して、フィティングカーブ81を生成し、ToFカーネルのスペクトルの同定を行った。続いて、得られたToFカーネルのスペクトルを基に、ML-EM法で再構成を行った。
図14に、再構成結果が示されている。グラフ90は、実施形態に係る方法による補正を行った場合の再構成結果、グラフ91は、実施形態に係る方法をによる補正を行わなかった場合の再構成結果を示している。グラフ90とグラフ91とを比較すると、ピーク付近でのカウント数が約2倍になっており、放射線物質の分布のピークが先鋭化することが確認できた。
【0172】
以上のように、実施形態では、発光モデルに基づいて、ガンマ線がシンチレータに入射してから発光が起こるまでの遅延時刻を推定し、推定結果に基づいてToFカーネルの補正を行う。このことで、ToFスペクトルを先鋭化することが可能となり、再構成後のPET画像の画質を向上させることができる。
【0173】
以上の実施形態に関し、発明の一側面および選択的な特徴として以下の付記を開示する。
【0174】
(付記1)
本発明の一つの側面において提供される核医学診断装置は、取得部と、算出部と、特定部と、測定部と、補正部とを備える。取得部は、第1の検出器で検出される第1の光子数情報を取得する。算出部は、前記第1の光子数情報に基づいて前記第1の検出器に対応する第1の発光確率モデルを算出する。特定部は、前記第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定する。測定部は、前記第1の検出器で検出されるイベントの検出タイミングを測定する。補正部は、前記第1のタイミングに基づいて、前記検出タイミングを補正する。
【0175】
(付記2)
前記第1の光子数情報は、シンチレーション光子数を含んでもよい。
【0176】
(付記3)
前記第1の光子数情報は、チェレンコフ光子数を更に含んでもよい。
【0177】
(付記4)
前記特定部は、前記第1の発光確率モデルにより示される分布の累積分布に基づいて、前記第1のタイミングを特定してもよい。
【0178】
(付記5)
前記特定部は、前記第1の発光確率モデルにより示される分布が極値を取る時刻に基づいて、前記第1のタイミングを特定してもよい。
【0179】
(付記6)
前記第1の発光確率モデルは、検出事象に含まれる複数の光子の内の最初の光子を検出する確率密度に関するものであってもよい。
【0180】
(付記7)
前記第1の発光確率モデルは、検出事象に含まれる複数の光子の内の所定の数の光子を検出する確率密度に関するものであってもよい。
【0181】
(付記8)
前記算出部は、前記第1の光子数情報に基づいて、各時刻における発光の分布関数を生成し、生成した前記分布関数に基づいて、前記第1の発光確率モデルを算出してもよい。
【0182】
(付記9)
前記取得部は、検出器毎に前記第1の光子数情報を取得してもよい。
【0183】
(付記10)
前記取得部は、イベントごとに前記第1の光子数情報を取得してもよい。
【0184】
(付記11)
前記算出部は、シンチレータ内での発光位置の不定性を考慮して、前記第1の発光確率モデルを算出してもよい。
【0185】
(付記12)
本発明の一つの側面において提供されるデータ処理方法は、第1の検出器で検出される第1の光子数情報を取得し、前記第1の光子数情報に基づいて前記第1の検出器に対応する第1の発光確率モデルを算出し、前記第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定し、前記第1の検出器で検出されるイベントの検出タイミングを測定し、前記第1のタイミングに基づいて、前記検出タイミングを補正する。
【0186】
(付記13)
本発明の一つの側面において提供されるプログラムは、第1の検出器で検出される第1の光子数情報を取得し、前記第1の光子数情報に基づいて前記第1の検出器に対応する第1の発光確率モデルを算出し、前記第1の発光確率モデルに基づいて、検出確率が所定の閾値以上となる第1のタイミングを特定し、前記第1の検出器で検出されるイベントの検出タイミングを測定し、前記第1のタイミングに基づいて、前記検出タイミングを補正する処理をコンピュータに実行させる。
【0187】
以上説明した少なくとも1つの実施形態によれば、画質を向上させることができる。
【0188】
いくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更、実施形態同士の組み合わせを行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【0189】
150a 取得機能
150b 算出機能
150c 特定機能
150d 測定機能
150e 補正機能
150f 再構成機能
150g 制御機能
150h 受付機能
150i 画像生成機能
150j 表示制御機能