JP6708561B2 - 放射線信号の処理方法、当該方法を実施する信号処理装置、コンピュータプログラム、当該プログラムを記憶したコンピュータ可読媒体 - Google Patents

Description

本発明は、信号を処理する方法、それに関連する信号処理装置、コンピュータプログラム要素、及びコンピュータ可読媒体に関する。

コンピュータトモグラフィは、通常、相対的に高い光子フラックスを有するＸ線ビームを必要とする。これは、このようなイメージングシステムにおけるエネルギー感受性光子計数検出器の使用を妨げることがあり、その理由は、高いフラックスが計数信号を歪めるパルスパイルアップ効果をもたらし、これが、そのような信号から再構成される画像にアーチファクトを生じさせるからである。しかしながら、光子計数検出器は、検出された信号の物質分解を実施することが望まれるスペクトルＣＴアプリケーションにおいて有用であることが分かっている。

今日の傾向は、一層高いフラックスを使用するスペクトルコンピュータトモグラフィスキャナの開発に向けられており、従って、パルスのパイルアップ効果は大きくなり、これは、実際、物質分解に関してＣＴアプリケーションを追及することをはばむ。

今日の文献において、多くのモデルが、パルスパイルアップを受ける期待されるカウント数を予測するために存在するが、この問題に対する汎用的な完全な解決法はまだ見つかっていない。例えばK Taguchi他の"An analytical model of the effects of pulse pileup on the energy spectrum recorded by energy resolved photon counting x-ray detectors", Med Phys 37(8), August 2010, pp 3957-3969を参照されたい。更に、パルスパイルアップ補正方法が記述される本願出願人の国際公開第2008/155679号公報を参照されたい。

識別される上述の問題に対処し、パルスパイルアップをより良く考慮するための方法又は装置のニーズがなおある。

本発明の目的は、独立請求項の本発明の主題によって解決され、他の実施形態は、従属請求項に組み込まれている。本発明の以下に記述される見地は、信号処理装置、コンピュータプログラム要素、及びコンピュータ可読媒体に等しく適用されることが留意されるべきである。

本発明の第１の見地によれば、放射線がサンプルを通過した後にエネルギー分解検出器システムによって測定される信号を処理する方法であって、検出器との光子相互作用に関して測定されるカウントレートを受け取るステップと、検出されたカウントに信号モデルをフィッティングするステップであって、前記サンプル中の物質の物理量が与えられる場合、モデルは、少なくとも１つのエネルギーレベル（Ｅ）の少なくとも１方向のパルス交差の頻度の条件付き期待値を表し、信号モデルが、パルス高さ変数（ｕ）のフーリエ表現に関して表現される、ステップと、を有する方法が提供される。

一実施形態によれば、信号モデルは、前記検出器システムの累積スペクトル応答関数（Φ_Ｒ）のフーリエ変換を組み込む。

一実施形態によれば、フーリエ変換は、事前に計算される。

一実施形態によれば、カウントレートは、１又は複数の電子パルスとして記録され、１又は複数の電子パルスは、それらがパルス高さを有する時点で検出され、前記信号モデルは、１又は複数のパルスに関するランダム過程の記述に基づき、ランダム過程モデルは、２つの確率密度から形成され、一方の確率密度は、２つの連続するパルス間の間隔の分布に関し、他方は、パルス高さの分布に関する。

一実施形態によれば、ランダム過程モデルは、ショットノイズのモデルである。

一実施形態によれば、フィッティング処理は、前記物理量のパラメータに関して前記モデルに基づいて目的関数を解くことを含む。

一実施形態によれば、物理量は、ｉ）吸収（Ａ）、ｉｉ）屈折、又はｉｉｉ）散乱の任意の１つである。

一実施形態によれば、物理量のパラメータに関して解く忠実性は、前記放射線の光子フルエンスから独立している。単なる統計エラー（ノイズ）はなお異なりうるが、「忠実性」は、バイアス又はシステマティックなエラーを意味する。

一実施形態によれば、信号モデルは、パルス形状の関数記述を含む。

第２の見地によれば、サンプルを放射線が通過した後に光子計数又はエネルギー分解検出器システムによって測定される信号を処理する方法であって、方法は、検出器との光子相互作用について測定されたカウントレートを受け取るステップと、検出されたカウントに信号モデルをフィッティングするステップであって、前記サンプル中の物質の物理量を与えられる場合、信号モデルが、少なくとも１つのエネルギーレベル（Ｅ）の少なくとも１方向のパルス交差の頻度の条件付き期待値を表し、フィッティング処理は、前記物理量のパラメータについて目標関数を解くことを含み、前記物理量のパラメータについて解かれる忠実性は、前記放射線の光子フルエンスから独立している、ステップとを有する方法が提供される。

一実施形態によれば、検出器システムは、光子計数タイプであり、又はエネルギー分解タイプであり、少なくとも１つのエネルギー閾値を規定することができ、レベル交差頻度は、少なくとも１つのエネルギー閾値の個々のものの出力カウントレートである。

要するに、我々がここに提案するものは、高いフラックスにおける光子計数検出器の完全な順方向モデルを構成する正確な分析式に基づく検出器信号処理方法である。

提案される方法は、基礎をなす信号モデルが、Ｘ線ビームに関して関心のある所与の物理量（例えば減衰、屈折又は出コヒーレンス／散乱能）について期待されるカウント数を予測するので、相対的に高いフラックスにおける光子計数検出器の効果的な動作を可能にする。特に、パイルアップ効果は、正しく考慮される。なお、「高いフラックス」に関して、典型的な臨床ＣＴ環境において、ｍｍ^２当たり毎秒１０^９のＸ線光子までのフラックス密度が、Ｘ線検出器のロケーションにおいて非減衰ビームで又はその付近で達成されることができる。

方法の例示的なアプリケーションにおいて、我々は、物質ベースの画像を導出するために、不正確なパイルアップ順方向モデルによってもたらされるバイアスに苦しまない最小２乗分解アプローチを提案する。

提案される方法は、特にスペクトルコンピュータトモグラフィに適用可能であり、この場合、全体の減衰を例えばヨウ素及び水又は骨及び軟組織のような異なる物質によってもたらされる減衰に分解するために、１又は複数（例えば５）のエネルギー閾値に測定データを量子化することが望まれる。提案される方法は、この分解を達成するために、正しい順方向又は逆方向モデを供給する。モデルベースの方法は、任意の検出器応答、任意に形作られた信号パルス及び任意のフラックスのために使用されることができ、更に正確な物質分解を可能にする。

本発明の例示的な実施形態は、添付の図面を参照して記述される。

スペクトルイメージング装置を示す図。 図１による装置の信号処理コンポーネントを示す図。 パルスパイルアップを示す図。 閾値ベースの光子計数を示す図。 信号処理方法のフローチャートを示す図。 シミュレーション結果と比較して、提案される方法による光子カウントの解析予測を示す図。 図６によるシミュレーションにおいて用いられるような放射線のフルエンススペクトルを示す図。

図１を参照して、スペクトルイメージング装置１００が示されている。スペクトルイメージング装置は、長手軸又はｚ軸を中心に検査領域２０８の周りを回転可能である回転可能ガントリ部分２０４を有するコンピュータトモグラフィ（ＣＴ）システムを有する。

Ｘ線管のようなＸ線源２１２は、回転ガントリ部分２０４によって支持され、検査領域２０８を横断するマルチエネルギーの放射線ビーム又は光子を放出する。

放射線感受性検出器２１６は、１又は複数のセンサ又はピクセル（「セル」）２１８を有し、ピクセルの各々は、検査領域２０８を横断する、線源２１２によって放出された光子を検出する。センサ２１８は、個々の検出された光子を示す電流又は電圧のような電気信号を生成する。適切なセンサの例は、例えば一般にシリコン、テルル化カドミウム（ＣｄＴｅ）又はテルル化カドミウム亜鉛（ＣＺＴ）から形成されるストリップのような半導体ウェハ部分又は本体を有する検出器のような直接変換を含む。他のオプションは、フォトセンサと光学通信するシンチレータを有するシンチレータベースのセンサである。

寝台のような対象支持体２４８は、検査領域２０８において患者又は他の対象を支持する。対象支持体２４８は、スキャニングプロシージャを実施する際に検査領域２０８に対して対象をガイドするように移動可能である。汎用コンピュータは、オペレーターコンソール２５２として有用である。コンソール２５２は、モニタ又はディスプレイのような人間可読の出力装置、及びキーボード及びマウスのような入力装置を有する。コンソール２５２に常駐するソフトウェアは、例えばグラフィカルユーザーインタフェース（ＧＵＩ）を通じてオペレータが、スキャナ２００を制御し、それと相互作用にすることを可能にする。このような相互作用は、より詳しく以下に説明されるように、エネルギービニングされたデータに基づいて信号を再構成する命令を有することができる。

イメージングセッション中、標本（人間又は動物の患者又は被検体の一部、又は他の関心対象（必ずしも有機でなくてよい）が、検査領域２０８に存在する。初期スペクトルを有するマルチエネルギーの放射線ビームが、標本を通過する。標本の通過中、放射線は、標本中の物質と相互作用し、その相互作用の結果として、変更される。この変更された放射線は患者から出て、検出器ピクセルと相互作用して対応する電気信号を生成する。多くの物理的なプロセスが、標本の物質を通るその通過中に、放射線の変化をもたらすように働く。その中でも特に、物理的なプロセスは、吸収、屈折、及び非コヒーレンス（小さい散乱角又は暗視野効果）である。それらの物理的プロセスの各々は、例えば局所吸収係数μ、局所屈折率φ、及び小角散乱能（small angle scattering power）Ωのような関連する物理量によって記述されることができる。前述の物理量の例えば吸収係数μは、それが概して標本全体にわたってその各々のポイントで異なるので、局所的である。特に、吸収ビューは、物質のタイプ（脂肪、骨又は他の物質）及びそのポイントにおけるその密度の関数である。更に、例えば吸収μにはエネルギー依存性がある。Alvarez及びMacovskiは、これに関して広範囲に記述している（"Energy-selective Reconstructions in X-ray Computerized Tomography", PHYS. MED. BIOL., 1976, VOL. 21, NO. 5, 733-744）。減衰又は吸収係数μが、物質の元素組成に特徴的である態様で、エネルギーによって変化することが知られている。言い換えると、放出されたＸ線エネルギースペクトルは、対象を通るその通過中、特徴的な「カラーリング（coloring）」を受ける。物質特有の減衰係数及び標本中の個々の物質を通る個々の線積分の線形結合として、経路に沿った、標本を通る全体の減衰が、表現されることができる。C Schirra他の"Statistical Reconstruction of Material Decomposed Data in Spectral CT", IEEE Trans on Medical Imaging, vol 32, No 7, July 2013のp1250のｅｑ（１）を参照されたい。これらのことが、物質ベースの関心のある要素の各々について特徴的な画像に至るためにスペクトルＣＴイメージングにおいて利用しようと試みられているものであり、全体のタスクは、検出信号を、さまざまな物質特有の線積分に分解することである。言い換えると、検出器において検出される電気信号は更に、詳しく後述する態様でスペクトル分解される。スペクトル処理されたデータは、再構成器２４４に伝送される。再構成器は、処理されたデータを、異なる物質の各々について特徴的な画像に再構成するように動作する。言い換えると、再構成器２４４が、センサ１１８によって生成された信号を、検出された光子のスペクトル特性に基づいて選択的に再構成するので、基本的な分解が達成される。例えば、ビニングされたデータは、概して、骨、有機組織、脂肪及び／又はその他のような異なる光子吸収特性を有する有機物質のタイプを分離するために、コントラスト強調物質を位置付けるために、及び／又は他のやり方でスペクトル特性に基づいて検出信号を処理するために、使用されることができる。標本は、Ｘ線源の回転の間、それぞれ異なる投影方向にまたがって放射線に曝露されるので、関心のある物質に関する内部の断面表現が再構成されることができる。

エネルギースペクトル情報は、パルス高さ解析ＰＨＡを通じて抽出される。ＰＨＡは、光子計数回路ＰＣＴを通じて図１の実施形態において実現される。ＰＣＴ回路は、検出器２１６においてさまざまな検出器セルと相互作用する。変更されたＸ線ビームの光子の各々は、特定のエネルギーを有する。このエネルギーは、電気パルスに正比例し、より具体的には、パルス高さに正比例しうる。より具体的には、特定の光子が、検出器セル２１８の（大きく空乏化された）半導体本体と相互作用する場合、ある数の電子正孔対が生成される。こうして生成された電子−正孔対の数は、相互作用する光子のエネルギーに正比例する。相対的に大きいバイアス電圧が、検出器セルの半導体本体の全体に印加される。バイアス電圧は、電子及び正孔が互いから離れて拡散されるようにする。正孔及び電子は、検出器セルの反対の側の個々の電極において収集される。この収集は、そのように拡散された電子−正孔対の数の関数として、ピクセル電極の読出しチャネルにおいて特徴的な電流又は電圧パルス高さを生じさせる。検出器セルにおいて生成される電子パルスは、以下の態様で光子計数回路ＰＣＴによって処理される。

プリアンプ２２０が、ピクセル２１８の任意のものによって生成される各々の電気信号を増幅する。

パルス整形器２２２は、検出された光子について増幅される電気信号を処理し、検出された光子を示す例えば電圧のようなパルス又は他のパルスを含む対応するアナログ信号を生成する。こうして生成されたパルスは、予め規定された形状又はプロファイルを有する。この例において、パルスは、検出された光子のエネルギーを示すピーク振幅を有する。

エネルギー弁別器２２４は、アナログパルスをエネルギー弁別する。この例において、エネルギー弁別器２２４は、アナログ信号の振幅を特定のエネルギーレベルに対応する閾値と比較する複数の比較器２２８を有する。言い換えると、弁別器２２４は、整形器２２２によって生成された入力パルスの「高さ」を決定するように動作する。より具体的には、各々の比較器２２８は、パルスの振幅がその閾値を越えるかどうかを示す出力信号を生成する。この例において、各々の比較器からの出力信号は、パルスの振幅が増大してその閾値を越えるとき、低から高（又は高から低）への遷移を含むデジタル信号を生成し、パルスの振幅が低下してその閾値を横切るとき、高から低（又は低から高）への遷移を含むデジタル信号を生成する。

例示的な比較器の実施形態において、各比較器の出力は、振幅が増大してその閾値を横切るとき、低から高に遷移し、パルスの振幅が低下してその閾値を横切るとき、高から低に遷移する。

カウンタ２３６は、閾値ごとにそれぞれに立ち上がりエッジ（又はある実施形態では立ち下がりエッジ）を計数する。カウンタ２３６は、閾値ごとに単一のカウンタ又は個別のサブカウンタを有することができる。エネルギービニング器２４０は、カウントを、エネルギー閾値の間のレンジに対応するエネルギーレンジ又はビンにエネルギービニングする。ビニングされたデータは、検出された光子をエネルギー分解するために使用される。

言い換えると、ＰＣＴ信号処理チェーンは、複数の電圧閾値によって規定されるエネルギービンに、各々の入力パルスのパルス高さを量子化するように動作する。Ｎ個の（Ｎは少なくとも１であり、好適には少なくとも２である）（電圧、電流量又は他のエネルギー標示）閾値が、規定されることができ、Ｎ個のエネルギービンが、前記閾値の個々のものより高いパルス高さを記録する。例えば、エッジが前記閾値のうち２つ（すなわち「十字」）を上回るパルスは、個々の２つの閾値に関連付けられる２つのビンの各々にカウントを生じさせる。閾値のより低いほうだけを横切る場合、ただ１つのカウントがある、などである。しかしながら、これは一例にすぎず、ある実施形態において、立ち下がりエッジのみがカウントを生じさせ、又は立ち上がりエッジ及び立ち下がりエッジの両方が、カウントを生じさせる。言い換えると、ＰＣＴ回路は、その出力において、各ビンにカウント数を供給し、いわばヒストグラムを与える。すなわち、ＰＣＴ回路は、変更されたＸ線ビームが検出器電子部品との相互作用において生じさせた、特定の時間フレームにおける全パルスの数及び高さを記録する。各ビン（１、２のエネルギービン又は好適には３又は５のようなより多くのエネルギービンが本願明細書では想定される）における記録された光子カウントは、スペクトルアナライザＳＰに伝送される。

スペクトルアナライザのコンポーネントが、図２に更に詳しく示される。概して、スペクトルアナライザＳＰは、入力ポートＩＮ及び出力ポートＯＵＴを有する。スペクトルアナライザＳＰは、２つのモジュール、予測器ＰＲ及びソルバ(solver)ＳＬＶを有する。予測器ＰＲ又はソルバＳＬＶによる高速フーリエ変換リクエストに答える高速フーリエ変換コンポーネントＦＦＴがある。フーリエ変換コンポーネントＦＦＴは、高速フーリエ逆変換を計算するために更に使用されることができる。更に、検出器電子部品の仕様(specifications)が記憶されるメモリＭＥＭがある。

予測器モジュールＰＲは、検出器仕様を読み取り、順方向信号モデルの関数記述にそれを含める。順方向信号モデルは、サンプル又は標本中の特定の物質の組み合わせが与えられる場合、各ビンごとの平均カウントレートを予測し、すなわち計算する。より詳しくは、サンプル中の１又は複数の異なる材料タイプの存在によってもたらされる関心のある物理的効果の１又は複数の処理が与えられる場合、予測器は、異なるエネルギー閾値の任意のものについて、前記閾値を横切るパルスを記録する条件付き期待値を計算する。言い換えると、条件付き期待値は、関心のある物理的効果の各々についてパラメータに関して表現される。前記パラメータは、フリーパラメータであり、光子カウントｃ＝（ｃ_１、ｃ_２、．．、ｃ_Ｎ）について実際に観察されたものに対するフィッティング処理において適応されることができる。より具体的には、観察された光子カウントは、入力ポートＩＮにおいて受け取られる。受け取られた光子カウントは、目的関数を公式化するためにモデルと組み合わせられる。目的関数は、ソルバＳＬＶによって、関心のある物理量パラメータについて解かれる。物理量Ａのパラメータについて解いたものが、出力ポートＯＵＴに出力される。

一実施形態において、こうして出力される物理的パラメータは、物質に特徴的な減衰係数を記述する。言い換えると、出力は、概して、所与の物質タイプの個々の物理量（例えば減衰係数）を記述する各エントリをもつベクトルである。概して、信号プロセッサＳＰは、所与の投影方向において取得された観察されたカウントレートから、関心のある物理量の前記方向に沿った投影へ、モデルベースの変換を実現する。ここでも、この物理量Ａは、物質特有の吸収係数でありえ、又は場合によっては屈折率又は散乱能でありうる。ここで提案される使用されるモデルは、単にパイルアップ挙動を近似するための文献における以前の試みを不要にする。ここで提案されるモデルは、任意の高いフラックスにおける光子計数検出器の完全な順方向モデルを構成する正確な解析演算を示す。言い換えると、ここで提案される解析モデルは、パイルアップ効果を正確にモデル化する。具体的には、パイルアップ効果がモデルの中で正しく考慮されるので、以前に使用されていたパイルアップ補正モジュールは必要ではない。当然ながら、これは、パイルアップ補正スキームにおいて、提案されたモデルベースのアプローチを使用することを除外せず、この形態は、他の理由で有益に使用されうる。

図３のグラフは、パイルアップ効果を示す。グラフは、検出器の半導体本体の光子相互作用に応じてパルス整形器によって出力される、それぞれのパルスのパルストレインを示す。明確さ及び単に例示の目的のために、同じ高さのパルスが示されているが、これは、当然ながら実際の例ではなく、なぜなら、パルス高さ自体が光子エネルギーによってランダムに変化するからである。パイルアップの出現は、円ＰＰに示される。パイルアップ事象は、約４００ナノ秒のところで生じた。検出器電子部品のパルス分解時間能力が許すよりも時間的に接近してパルスが到達する場合に、パイルアップが生じる。言い換えると、図から分かるように、第２の、遅いほうのパルスは、早いほうのパルスのテイルに「乗り」、それにより、組み合わされたパルス（合成パルス）を形成する。しかしながら、合成パルスは、各パルスから組み合わされたものであり、各パルスよりも高くなっている。これは、「二重の」偽読み取りをもたらす：ｉ）不正確な数のパルスが記録される（２つのパルスが一体となって記録される）；ｉｉ）不正確なパルス高さが記録される（パルス高さが、「誤った」重ね合わせのため過大に評価される）。これは、このようなパルスに影響されるカウントデータ（更に本明細書において、ＯＣＲ（すなわち「出力カウント−レート」）と称される）から再構成される画像に、深刻なアーチファクトをもたらすことがある。

図４は、より現実的な例を示しており、すなわち、パルスの各々の出現は、光子相互作用事象に対応しており、各パルスのパルス高さは、光子のエネルギーによって異なる。更に、図４は、カウントｃがビニングされる３つの異なるエネルギービンを規定するためにＤＣ、ＤＣの０．５倍、及びＤＣの１．５倍と示される３つの閾値レベルに関して、入力パルスをインタロゲートするＰＣＴ電子部品の事例を示す（０．５以下のパルスは検出されない）。図から分かるように、記録されるカウントは、１００Ｍｈｚの高いレートのため大きく歪められ、従って、多くのパルスがパイルアップする。ここでも、任意の数の１又は複数の閾値が企図され、図４では３つの閾値が示されている。

図３及び図４のグラフは、パルスの時間連続又は時間トレインを示す。しかしながら、時間連続自体は、検出器によって観察されず、観察されるのはカウントイベントのみである。それでもやはり、ここに示すように、パルスの「ピーク」を、例えば図３、図４の電圧−時間のダイアグラムにマップされる要素としてとらえ、それらのピークポイントの「動態」を考えることが有用であることが分かる。ピークポイントの動態は、二重確率プロセスＸ（ｔ）としてモデル化されることができる。プロセスは、それが、ｉ）時間軸に沿ったパルス間間隔（すなわち、２つの連続するパルスイベントの間の間隔）の分布をモデル化する１つの確率密度と、ｉｉ）前記間隔から独立して、それ自身が独立確率変数であるパルス高さＵをモデル化する第２の確率密度と、を使用するので、「二重」確率である。電圧−時間面におけるピークポイントの動態を考えることが可能である。そして、「速度」「Ｘドット」をピークポイントに起因するものとすることができる。さらに、単位時間あたりの前記ピークポイントのレベル交差を検討することが可能である。更に発見的手法を利用することにより、各レベルの周辺に近傍を規定することができる。前記近傍内のピークポイントの速度が分かっている場合、速度を、個々のレベルの期待される交差と関連付けることができる。これと同様の式は、ショットノイズを研究したO. Riceによって１９４０年にすでに提案されており、O. Riceは、特に、動的過程×（ｔ）及びその時間微分ｔである「Ｘドット」の以下の同時確率分布を含む静止エルゴード過程のレベル交差頻度を研究した：

S. O. Rice, Mathematical Analysis of Random Noise, Bell Systems Tech. J., Volume 24, p. 46-156, 24:46 -156, 1945を参照されたい。

本願出願人は、ショットノイズの確率動力学が、光子計数装置の検出器電子部品のカウントイベントを正確に記述することを見つけた。本願出願人は、パイルアップ現象がショットノイズと数学的に同一であり、エネルギー感受性Ｘ線検出器の計数結果が、ショットノイズプロセスのレベル交差頻度によって識別されることができることを実験によって見出した。より最近の数学の刊行物（H Bierme and A Desolneux, "A Fourier Approach for the Level Crossings of Shot Noise Processes with Jumps", Journal of Applied Probability, 49(1):100 -113, 2012）は、ショットノイズプロセスのレベル交差が、頻度自体を考えるのではなく、レベル交差頻度のフーリエ変換を考える場合に、数学的により扱いやすくなることを示している。フーリエ空間でのこの簡略化を利用することにより、パルス高さＵの空間へのフーリエ逆変換（最も外側の積分を参照）後の順方向（又は逆方向）信号モデルの解析表現は、以下のように書くことができる：

上式において、Ｎ_ｘ（Ｕ，Ａ）は、所与のエネルギー閾値Ｕ及び対象による所与の減衰Ａについて記録されたカウントの期待されるレートを示す。Ｄ（Ｅ）は、検出器の吸収効率であり、Φ_０（Ｅ）は、Ｘ線源スペクトルであり、

は、累積検出器応答関数のフーリエ変換を示し、μ_α（Ｅ）は、基礎物質αの線減弱係数であり、「パイルアップ」カーネル、

は、以下によって与えられる：

上述のパイルアップカーネル（２）は、エネルギーＥにおける任意の数のパルスの複雑な線形の重ね合わせを正確に記述しており、検出器信号のアナログ出力が単位時間当たりのパルス高さＵについて閾値を超えた回数及びパイルアップによる変更を予測する。

カーネル演算（２）が行うことは、累積検出器応答のフーリエ変換によって各々加重されるパルス波形を、時間にわたって積分することである。式（１）において、大文字のＵは、パルス高さ空間（例えば電圧を表す）を示し、小文字ｕは、以下のフーリエ規定のようにフーリエ空間においてパルス高さＵの頻度変数を示す：

上式で、Ｎ_Ｘ（Ｕ）は、プロセスＸ（ｔ）によるレベルＵの期待される交差レートを示す。言い換えると、式（１）のモデルは、フーリエ空間において変数ｕに関して表現される。ｕは、Ｕに関連するフーリエ空間における頻度変数を表す。

パルスピークのシーケンスの基礎をなす二重確率モデリングは、以下によって与えられる：

上式で、Ｘは、すでに示された二重ランダム過程の確率変数であり、ｔ_ｉは、時間的なランダム指数ポアソン分布に続くポアソン点であり、パルス高さ確率変数Ｕは、パルス高さの空間においてスペクトルｐ（Ｕ）を有し、Ｕ_ｉは、異なるｉについて独立しているが等しく分布する。量ｇは、パルス形状を「インプリント（imprint、特徴づけ）」し、ＰＣＴのパルス整形器の電子部品によって決定される。例示のパルス形状は、g(t)=e^-t/τ(t>0)であり、τは、典型的なパルス持続時間であるが、長さａの矩形形状のような他のパルス波形も可能であり、高さＵ_０が更に可能である。

Bierm他の結果を拡張することにより、本願出願人は、パイルアップによるパルスの複雑な重ね合わせが、パルス高さスペクトルｐ（Ｕ）のフーリエ変換、

を考えることによって、フーリエ空間に「巻き込まれることができる」ことを見出した。Ｎ_Ｘ（Ｕ）のフーリエ変換は、立ち上がりパルスエッジが生じるとき（ジャンプするとき）の、フーリエ変換されたスペクトルｐ（ｕ）の応答の合計として表現されることができることが分かる。前記ジャンプの際の応答の合計は、同時確率密度、

のフーリエ変換を組み込む速度成分によって重み付けられる。Biermの第６頁の式７を参照されたい。

光子計数の目的での検出器電子部品の挙動に対するショットノイズコンテクストの適用は、パルス高さ密度ｐ（Ｕ）が、

のように、検出器応答Ｒ（Ｅ，Ｕ）を通じて表現されることができるという本願出願人の認識の部分において導出されている。

これを、上述の

のようにフーリエ空間に変換することにより、パルスパイルアップの重ね合わせの望ましい簡略化を与える。

上述された形の式（１）は、記録されるカウントの期待されるレートの減衰ベース｛μ_α（Ｅ）｝の依存及び対応する基礎物質の線積分Ａ＝（Ａ_α）を明らかにする。これは、mAs sr keV当たりの光子の単位でソースの光子フルエンスΦ_０を導入することによって達成され、ここで、ｓｒは、ステラジアンである。言い換えると、式１に従う信号順方向モデルは、関心のある物理量、この例では減衰について、期待されるＯＣＲを条件づけることによって、書き出された。同様の式は、関心のある他の物理量の各々についても、すなわち位相コントラストイメージングの場合の屈折φ及び暗視野イメージングの場合のデコヒーレンス散乱についても、導出されることができる。更に、パルス高さは、ここでは電圧として表現されているが、これが他の状況において必ずしもそうでなくてもよく、代わって他の適切な観測可能な電流量が、ＯＣＲを表現する又は定量化するために使用されることができる。

式１の構造から、それが２つのフーリエ積分の積のフーリエ逆変換（これは外側の積分）を含むことを理解することができる。信号モデルの各々の成分は、累積応答関数Φ_Ｒ（Ｅ，Ｕ）のフーリエ変換、

を含み、累積応答関数Φ_Ｒ（Ｅ，Ｕ）は、Ｕより大きいパルス高さを有するエネルギーＥのイベントごとに記録されるパルスの平均数を評価する。累積応答関数Φ_Ｒ（Ｅ，Ｕ）は、Schlomka他の"Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography", Phys. Med. Biol. 53 (2008) 4031 -4047 on page 4037 in section 4.1によって記述されるように、単色入力Ｘ線の低いフルエンスにおいて、閾値スキャンの間に容易に測定されることができる。Schlomkaの文献は、参照によってその全体がここに盛り込まれるものとする。

式１による信号モデルの構造は、主な入力コンポーネントが、検出器電子部品の仕様であることを示す。より具体的には、検出器吸収効率Ｄ及びＸ線源スペクトルΦが知られると、特定することが残っているものは、累積応答関数（のフーリエ変換）、

及び図１のより早いパルス整形器のようにＰＣＴ回路のパルス整形器によって出力されるパルス形状関数ｇ（ｔ）だけである。これらの量が与えられると、モデル（１）は、関心のある物理量についてセットアップされることができ、例えばこの例では、物理量は、減衰であり、個々の物理パラメータは、各基礎物質αの線積分Ａである吸収係数の投影である。

式（１）において極限をとる解析は、低い光子フルエンス（又は高減衰、又はそれらの両方）について、

を含む指数項が１に近づき、低フラックス期待値が、

のように回復されることを表す。

式（１）の構造は、基本的に、２つの関数の積のフーリエ変換の形であり、従って、式（１）は、「低フラックス」期待値の式（３）及び式（１）の指数項のフーリエ変換（後者は、低フルエンスに関するデルタ関数になる）のパルス高さ空間における畳込み、として書き直されることができ、言い換えると、

ここで、

信号プロセッサＳＰの動作が、図５のフローチャートを参照して説明される。

ステップＳ５０５において、検出器電子部品の仕様が読み込まれる。検出器の仕様は、特に累積応答関数と、光子計数データｃを取得するために使用される検出器電子部品のパルス整形器によって供給されるパルス形状の関数記述と、を含む。累積応答関数及び／又はパルス形状記述はフーリエ変換される。フーリエ変換された検出器の仕様は、現在使用されている検出器電子部品に関して初期化するために、上述の式１のモデルに含められる。

モデルは、物質分解処理において再構成することが望まれる物理量を反映する関心のある適切な物質特有のパラメータを含むように、順方向モデルとしてセットアップされることができる。一実施形態により、物質特有の減衰係数から形成されるパラメータＡは、減衰係数を特定し、又はサンプルを通る個々のＸ線源−検出器ピクセル経路に沿った個々の線積分が、さまざまな基礎物質から作られることが考えられる。

ステップＳ５１０において、実際のＯＣＲカウント（ビン当たりのカウント）が受け取られる。

ステップＳ５１５において、フィッティング処理が、モデルの物理量パラメータを、受け取ったＯＣＲカウントにフィットさせるために実施される。一実施形態において、図１のモデル（１）の予測されるレベル交差レートが、目的関数を通じて、実際に観察されたカウントレートｃに関連付けられる。式（１）のスペクトル順方向モデルの１つの実現例において、フーリエ変換表現、

が、事前に計算されることができ、それらは関心のある減衰又は個々の物理量から独立しているので、要求時に取り出すためにメモリＭＥＭに保持されることができる。特定の減衰Ａが与えられる場合に、（各閾値Ｕについて）期待されるカウント数を計算するために、（フーリエ空間のあらゆるｕについて）エネルギーに関する２つの線積分及び１つのフーリエ変換が必要とされる。しかしながら、式（１）のエネルギーに関する第１の積分は、フーリエ変換の後まで延期されることができ、これは、計算時間を節約する。更に、フーリエ変換は、Ｕの値ごとに計算される必要はなく、所与の状況において関心のある閾値についてのみ必要とされる。言い換えると、提案される方法は、関心のある閾値（すなわち例えばＰＣＴ電子部品において実際にセットアップされた閾値）を確立するための事前選択ステップを含むことができ、従ってＦＦＴコンポーネントが指示される。

フィッティングステップＳ５１５の出力は、受け取ったビンカウントから、所与の投影方向及び／又は経路についてセットされた物理パラメータへの変換を一緒に規定する物理パラメータの組である。

フィッティング処理Ｓ５１５は、投影方向及び検出器素子ごとに実施される。ステップＳ５１５の目的関数は、例えば最小二乗法乗又は最尤法（例えばガウス又はポアソンのような適切なノイズモデルを必要とする）等の任意の適切な最適化アルゴリズムによって、所望の物理パラメータについて解かれることができる。

より具体的には、フィッティングステップＳ５１５は、放射線源から個々の検出器ピクセルまでの各々の可能性のある経路について繰り返される。言い換えると、基礎物質ごとの物理パラメータの組が、各検出器ピクセルについて計算される。

ピクセルと経路との間の関連付けは、イメージャの制御部品によって供給され、イメージャの制御部品は、変化するイメージングジオメトリを追跡し、カウントの各々の組を個々の投影経路で特徴付ける。こうして計算された物理パラメータＡは、ステップＳ５２０に伝送され、ステップＳ５２０において、フィルタード逆投影その他の再構成アルゴリズムが、基礎物質ごとに物理パラメータの各々に適用され、複数の異なる断面画像が得られ、各物質について１の断面画像がある。以下に、物理量が基礎物質ごとの吸収である例が与えられる。

例示的な実施形態において、我々は、以下のように上述の順方向モデルに基づく分解技法を提案する。ｍ_ｉを、閾値Ｕ_ｉについての測定されたカウント値とし、Ｕ_１...Ｕ_Ｍを、エネルギー感受性光子計数検出器のＭ個の異なるエネルギー閾値とする。次に、我々は、Ａに関して以下の関数（又は目的関数）を最小にすることによって、以下の最小２乗アプローチから物質の線積分を評価することができる：

最小２乗推定は、以下によって与えられる：

上式において、σ^２ _ｉ（Ｕ_ｉ，Ａ）は、測定値ｍ_ｉの分散であり、それらは、期待値自体がそうであるように、概してパルス高さ閾値Ｕ_ｉ及び減衰Ａに依存する。しかしながら、この効果は二次的であり、分散についての簡略化されたモデルは、物質の線積分を信頼性をもって評価するのに十分でありえ、例えばポワソン近似は次のようになる：

例えばガウス又はポワソンノイズモデルを用いる最尤法アプローチのような、最小二乗アプローチに代わる方法が、当然ながら潜在的な推定器として可能である。

図１による信号モデルは、ここでは順方向モデルとして使用され、言い換えると、カウントは、関心のある物理量（吸収、屈折又は暗視野散乱）が与えられる場合に予測されることが分かるであろう。式１のモデルは、カウントが与えられる場合に物理量が予測される逆方向モデルとしても見ることができる。更に、式１は、関心のある特定の物理量、すなわち吸収、に依存して書かれている。しかしながら、これは単に例示のためであり、例えば位相コントラストイメージングの場合の屈折及び暗視野イメージング場合のデコヒーレンスのような関心のある物理量に関して、式１のモデルを書き直すことは表記の問題である。提案されるモデルの正確さは、実質的に任意のフルエンスでのカウント、パルス高さの解析を可能にする。より具体的には、上述の方法を使用して所与の物理量について基準画像を計算し、この画像を異なるフルエンスレートで処理された他の画像のシーケンスと比較する場合、それらの画像についての平均バイアスは、両方のケースにおいて最小限にされることが分かるであろう。これは、期待されるカウント数に関しての単なる近似であるモデルが使用される場合には当てはまらない。これは、提案されたモデリングスキーム及びフルエンススペクトル全体にわたる信号処理方法の有用性を実際に示す。

提案された方法の正確さを、Matlab suite計算を通じて得られる数値データによって補強するために、図６が示される。より具体的には、モンテカルロ（ＭＣ）シミュレーションが、異なるデータセットを生成するために実行された。理解されるように、上述のモデルを使用した解析予測及びシミュレーションの結果は、実質的に一致する。曲線は、垂直軸上のレベル交差頻度レート（Ｍｃｐｓ）対水平軸上の電流（ｍＡ）を示す。ディスクマーカは、シミュレーションの進行を示し、実線は解析結果を示す。より具合的には、エネルギー閾値Ｕ＝３０ｋｅＶ、５０ｋｅＶ、７０ｋｅＶ、９０ｋｅＶ及び１１０ｋｅＶの５つの異なる値についての、光子計数検出器において記録される期待されるカウント数のモデルベースの解析予測及びモンテカルロシミュレーションが、Ω＝１０^−６ｓｒの立体角をスパンする検出器について、アノード管電流の関数としてプロットされる。

エラーバーは、Schlomkaの文献のＦＷＨＭによるガウス検出器応答（すなわち、ＦＷＨＭ＝ａ＋ｂＥ、ここでａ＝１：６１５ｋｅＶ、ｂ＝０：０２５であり、Ｅは、ｋｅＶで表されるエネルギーである）について、ＭＣシミュレーションの統計誤差（シミュレーションにおいて１ｍＶが１ｋｅＶに対応する）の５つの標準偏差を表す。Ｘ線フルエンススペクトルが、図７に示される。パルス波形は、上述したように指数形式であり、正のｔの場合はｇ（ｔ）＝ｅ^−ｔ／τであり、他の場合はゼロであり、ここで、τ＝１０ｎｓである。

ノイズ実現の数及び実現及びデータポイント当たりのパルスの数は共に１０００に選ばれた。エラーバーは、１０００の実現にわたるサンプル分散から得られるＭＣシミュレーションの５つの標準偏差を表現する。残りの偏差は、非常に小さく、有限の長さに関する離散的フーリエ変換の使用に起因する。

要約すると、我々がここで提案するものは、任意の検出器応答、パルス波形及び高いＸ線フルエンスに関するエネルギー感受性光子計数検出器の解析的な順方向又は逆方向モデルである。単位時間あたりの期待されるカウント数は、フーリエ積分によって、又は、低フルエンス期待及びパイルアップカーネルの畳込みによって記述される。図６に関するモンテカルロシミュレーションは、我々の結果の有効性を検証する。

上述の方法ステップは、物理量が物質特有の減衰係数である例において説明されたが、これは単に説明のためにすぎず、方法は、スペクトル位相コントラスト又はスペクトル暗視野イメージングについて容易に適応されることができる。例えば、位相シフト及び／又は散乱能に関してモデル（１）を考慮するために、吸収モデル項、

にモデル式、

を乗じることによって、モデルの複雑さが拡張される。
ここで、Ｖ（Ｅ）は、（干渉計の）可視性であり、Ωは、見通し線に沿った散乱能であり、ｆ_ＤＣは、のちの小角散乱のエネルギー依存性であり、Δφは、屈折によってもたらされる位相シフトに関し、αは、エネルギー依存のファクタである。Ｆ（）は、位相ステッピングにおいて又は任意の他の干渉計方法においてサンプル誘導される屈折によってもたらされる輝度変調をモデル化するための三角関数（例えばシヌソイド）であり、位相ステッピングは、１つの例示の実施形態である。式（１）のエネルギー積分は、拡張されたモデルのエネルギー依存項の各々を通じて拡張される必要がある。

式（１）のモデルにおいて、単極パルス構造が仮定されている。言い換えると、パルスは、常に正であり又は常に負である。しかしながら、式（１）は、正及び負の部分の両方をもつ双極パルスを考慮するために改良されることもできる。概して、上述のアプローチは、ステップ関数によって近似されることができるパルスに有効である。

更に、上述のモデルにおいて、カウントは、立ち上がりパルスエッジのみを記録することによって規定されるが、モデルは、立ち下がりエッジについても同等に適用される。

本発明の別の例示的な実施形態において、上述の実施形態のうちの１つによる方法の方法ステップを適当なシステム上で実行するように適応されることによって特徴付けられるコンピュータプログラム又はコンピュータプログラム要素が提供される。

従って、コンピュータプログラム要素は、コンピュータユニットに記憶されることができ、かかるコンピュータユニットもまた、本発明の実施形態の部分でありうる。このコンピューティングユニットは、上述の方法における各ステップを実施するように又はその実施を引き起こすように適応されることができる。更に、コンピューティングユニットは、上述の装置のそれぞれのコンポーネントを動作させるように適応されることができる。コンピューティングユニットは、自動的に作動するように及び／又はユーザの命令を実行するように、適応されることができる。コンピュータプログラムは、データプロセッサの作業メモリにロードされることができる。こうして、データプロセッサは、本発明の方法を実施する能力をもつことができる。

本発明のこの例示的な実施形態は、最初から本発明を使用するコンピュータプログラム、及び更新によって既存のプログラムを本発明を使用するプログラムにするコンピュータプログラム、の両方をカバーする。

更に、コンピュータプログラム要素は、上述した方法の例示的な実施形態のプロシージャを実現するためのすべての必要なステップを提供することが可能である。

本発明の他の例示的な実施形態により、コンピュータ可読媒体（例えばＣＤ−ＲＯＭ）が提示され、かかるコンピュータ可読媒体は、それに記憶されるコンピュータプログラム要素であって、上述のセクションによって記述されたコンピュータプログラム要素を有する。

コンピュータプログラムは、例えば他のハードウェアと共に又はその一部として供給される光学記憶媒体又はソリッドステート媒体のような適切な媒体に記憶されることができ及び／又は配布されることができるが、他の形態において、例えばインターネット又は他のワイヤード又はワイヤレス通信システムを通じて、配布されることもできる。

しかしながら、コンピュータプログラムは、ワールドワイドウェブのようにネットワークを通じて提供されることもでき、このようなネットワークからデータプロセッサの作業メモリにダウンロードされることができる。本発明の他の例示的な実施形態により、本発明の上述の実施形態のうちの１つによる方法を実施するように適応されるコンピュータプログラム要素をダウンロード可能にする媒体が提供される。

本発明の実施形態は、それぞれ異なる発明の主題に関して記述されていることに留意されたい。特に、ある実施形態は、方法タイプの請求項に関して記述され、他の実施形態は、装置タイプの請求項に関して記述されている。しかしながら、当業者であれば、上述及び後述の説明から、他の場合が示されない限り、１つのタイプの発明の主題に属する特徴同士の任意の組み合わせに加えて、個々の異なる主題に関連する特徴同士の任意の組み合わせが、本願によって開示されていると考えられることが分かるであろう。しかしながら、すべての特徴は、特徴の単純な足し合わせより多くの相乗効果を提供して組み合わせられることができる。

本発明は、図面及び上述の記述において詳細に図示され記述されているが、このような図示及び記述は、制限定なものではなく、説明的又は例示的であると考えられることができる。本発明は、開示された実施形態に制限されない。図面、開示及び従属請求項の検討により、開示された実施形態に対する他の変更が、請求項に記載の本発明を実施する際、当業者によって理解され達成されることができる。

請求項において、「含む、有する（comprising）」という語は、他の構成要素又はステップを除外せず、不定冠詞「a」又は「an」は、複数性を除外しない。単一のプロセッサ又は他のユニットは、請求項に列挙されるいくつかのアイテムの機能を果たすことができる。特定の手段が相互に異なる従属請求項に列挙されているという単なる事実は、これらの手段の組み合わせが有利に使用されることができないことを示さない。請求項における参照符号は、その範囲を制限するものとして解釈されるべきでない。

Claims (14)

1. 放射線がサンプルを通過した後にエネルギー分解検出器システムによって測定される信号を処理する方法であって、
   前記エネルギー分解検出器システムとの光子相互作用について測定されたカウントレートを受け取るステップと、
   前記検出されたカウントレートに信号モデルをフィッティングするステップと、
   を有し、
   前記信号モデルは、パルスパイルアップ効果を考慮するように構成され、前記信号モデルが、前記サンプル中の物質の物理量が与えられる場合に、少なくとも１つのエネルギーレベルの少なくとも一方向のパルス交差の頻度の条件付き期待値を表し、
   前記信号モデルは、パルス高さ変数のフーリエ変換を組み込む、方法。
2. 前記フィッティングの処理は、前記物理量のパラメータについて前記信号モデルに基づく目的関数を解くことを含む、請求項１に記載の方法。
3. 前記信号モデルは、前記エネルギー分解検出器システムの累積スペクトル応答関数のフーリエ変換を組み込む、請求項１又は２に記載の方法。
4. 前記カウントレートは、１又は複数の電子パルスがパルス高さを有する時点で検出される該１又は複数のパルスとして記録され、
   前記信号モデルは、前記１又は複数のパルスのランダム過程の記述に基づき、
   前記ランダム過程に基づく前記信号モデルは、２つの連続するパルス間の間隔の分布に関する確率密度と前記パルス高さの分布に関する確率密度との２つの確率密度から形成される、請求項１乃至３のいずれか１項に記載の方法。
5. 前記物理量が、ｉ）吸収係数、ｉｉ）屈折率、及びｉｉｉ）散乱能のうちの１つである、請求項１乃至４のいずれか１項に記載の方法。
6. 前記フーリエ変換が事前に計算される、請求項３に記載の方法。
7. 前記ランダム過程に基づく前記信号モデルが、ショットノイズのモデルである、請求項４に記載の方法。
8. 前記物理量のパラメータについて解くことの忠実性は、前記放射線の光子フルエンスから独立している、請求項１乃至７のいずれか１項に記載の方法。
9. 前記信号モデルが、パルス形状の関数記述を含む、請求項１乃至８のいずれか１項に記載の方法。
10. 放射線がサンプルを通過した後に光子計数検出器システムによって測定される信号を処理する方法であって、
    前記光子計数検出器システムとの光子相互作用について測定されたカウントレートを受け取るステップと、
    前記検出されたカウントレートに信号モデルをフィッティングするステップと、
    を有し、
    前記信号モデルは、前記サンプル中の物質の物理量が与えられる場合、少なくとも１つのエネルギーレベルの少なくとも一方向のレベル交差の頻度の条件付き期待値を表し、
    前記フィッティングの処理は、前記物理量のパラメータについて目的関数を解くことを含み、
    前記物理量のパラメータについて解くことの忠実性は、前記放射線の光子フルエンスから独立している、方法。
11. 前記光子計数検出器システムは、光子計数タイプであり、少なくとも１つのエネルギー閾値を規定することができ、
    前記レベル交差の頻度は、少なくとも１つのエネルギー閾値の個々についての出力カウントレートである、請求項１０に記載の方法。
12. 請求項１乃至９、１０又は１１のいずれか１項に記載の方法を実施する信号処理装置。
13. データ処理ユニットによって実行されるとき、請求項１乃至９、１０又は１１のいずれか１項に記載の方法の各ステップを実施するように適応される信号処理装置を制御するためのコンピュータプログラム。
14. 請求項１３に記載のコンピュータプログラムを記憶したコンピュータ可読媒体。

Images