以下に、図面を参照しながら、本発明の実施形態に係るX線CT装置、画像処理装置およびプログラムを詳細に説明する。また、以下の図面において、同一の部分には同一の符号が付してある。ただし、図面は模式的なものであるため、具体的な構成は以下の説明を参酌して判断すべきものである。
(第1の実施形態)
図1は、第1の実施形態に係るX線検査装置の全体構成図である。図1を参照しながら、X線検査装置1の全体構成の概要を説明する。
X線CT装置の一例であるX線検査装置1は、図1に示すように、放射線の一例であるX線を被検体40に透過させてエネルギーごとのフォトン数で示されるスペクトルとして検出することにより、被検体40の投影断面41の断面画像を得るスペクトラルCT装置またはフォトンカウンティングCT装置等である。X線検査装置1は、図1に示すように、架台装置10と、寝台装置20と、コンソール装置30(画像処理装置)と、を備えている。
架台装置10は、被検体40に対してX線を照射して透過させ、上述のスペクトルを検出する装置である。架台装置10は、X線管11と、回転フレーム12と、検出器13と、照射制御部14と、架台駆動部15と、データ収集部16(収集部)と、を備えている。
X線管11は、照射制御部14から供給される高電圧によりX線を発生する真空管であり、被検体40に対してX線ビーム11aを照射する。X線管11から照射されるX線のエネルギーごとのフォトン数で示されるスペクトルは、X線管11の管電圧、管電流、および、線源に用いるターゲット(例えば、タングステン等)の種類によって定まる。そして、X線管11から照射されたX線は、被検体40を透過する際に、被検体40を構成する物質の状態に応じてX線のエネルギーが減弱し、各エネルギーのフォトン数が減少してスペクトルが変化する。
回転フレーム12は、X線管11と検出器13とを、被検体40を挟んで対向するように支持するリング状の支持部材である。
検出器13は、チャネル毎に、X線管11から照射されて被検体40を透過したX線であるX線ビーム11bのエネルギーごとのフォトン数を検出する検出器である。すなわち、検出器13は、チャネル毎に、後述する図4に示すようなX線のエネルギーごとのフォトン数で示されるスペクトルを検出する。ここで、検出器13により検出されたスペクトルを、以下、「検出スペクトル」という場合があるものとする。検出器13は、図1に示すように、回転フレーム12の周方向に回転しながら、ビュー毎にスペクトルを検出する。ここで、ビューとは、回転フレーム12の周方向の1周360°のうち、所定角度ごとに検出器13によりスペクトルが検出される場合の角度のことをいうものとする。すなわち、検出器13が0.5°ごとにスペクトルを検出する場合、1ビュー=0.5°というものとする。検出器13は、チャネル方向(回転フレーム12の周方向)に複数の検出素子が配列された検出素子列が、被検体40の体軸方向(スライス方向)(図1に示すZ軸方向)に沿って複数列配列された2次元アレイ型検出器である。なお、検出器13の検出素子列は、フォトカウンティング型検出素子と積分型検出素子との組み合わせで構成されていてもよい。また、X線管11と検出器13との組が、複数組設置されているものとしてもよい。
照射制御部14は、高電圧を発生して、発生した高電圧をX線管11に供給する装置である。
架台駆動部15は、回転フレーム12を回転駆動させることで、被検体40を中心とした円軌道上でX線管11および検出器13を回転駆動させる装置である。なお、架台駆動部15は、X線管11および検出器13の双方を回転駆動させる構成に限定されるものではない。例えば、検出器13は、回転フレーム12の周方向に1周分にわたって検出素子が配列されて構成されており、架台駆動部15は、X線管11のみを回転駆動させる構成であってもよい。
データ収集部16は、検出器13によりチャネル毎に検出されたエネルギーごとのフォトン数で示されるスペクトル(第1スペクトル)のデータを収集する装置である。そして、データ収集部16は、収集したスペクトルのデータそれぞれに対して増幅処理またはA/D変換処理等を行なって、所定幅のエネルギー帯ごと(以下、単に「エネルギーごと」という場合もあるものとする)のサイノグラムを生成し、コンソール装置30に出力する。
寝台装置20は、被検体40を載せる装置であり、図1に示すように、寝台駆動装置21と、天板22とを、備えている。
天板22は、被検体40が載置されるベッド等の寝台である。寝台駆動装置21は、天板22に載置される被検体40の体軸方向(Z軸方向)へ移動させることによって、被検体40を回転フレーム12内に移動させる装置である。
コンソール装置30は、操作者によるX線検査装置1に対する操作を受け付け、架台装置10によって収集されたデータから断面画像(復元画像)を再構成する装置である。コンソール装置30は、図1に示すように、入力装置31(入力部)と、表示装置32と、スキャン制御部33と、画像処理部34と、画像記憶部35と、システム制御部36と、を備えている。
入力装置31は、X線検査装置1を操作する操作者が各種指示を操作入力するための装置であり、操作入力された各種コマンドをシステム制御部36に送信する装置である。入力装置31は、例えば、マウス、キーボード、ボタン、トラックボール、またはジョイスティック等である。
表示装置32は、入力装置31を介して操作者から操作指示を受け付けるためのGUI(Graphical User Interface)を表示したり、後述する画像記憶部35が記憶する復元画像(断面画像)を表示する装置である。表示装置32は、例えば、CRT(Cathode Ray Tube)ディスプレイ、LCD(Liquid Crystal Display:液晶ディスプレイ)、または有機EL(Organic Electro−Luminescence)ディスプレイ等である。
スキャン制御部33は、照射制御部14、架台駆動部15、データ収集部16、および寝台駆動装置21の動作を制御する処理部である。具体的には、スキャン制御部33は、回転フレーム12を回転させながら、X線管11からX線を連続的または間欠的に出射させることで、X線スキャンを実行させる。例えば、スキャン制御部33は、天板22を移動させながら回転フレーム12を連続回転させて撮像を行なうヘリカルスキャン、または、被検体40の周りを回転フレーム12が1回転して撮像を行い、続いて被検体40が載置された天板22を少しずらして再び回転フレーム12が1回転して撮影を行うノンヘリカルスキャンを実行させる。
画像処理部34は、データ収集部16から受信したサイノグラムから被検体の断面画像を再構成する処理部である。画像処理部34のブロック構成および動作の詳細については、後述する。
画像記憶部35は、画像処理部34による再構成処理により生成された断面画像(復元画像)を記憶する機能部である。画像記憶部35は、例えば、HDD(Hard Disk Drive)、SSD(Solid State Drive)、または光ディスク等の記憶装置である。
システム制御部36は、架台装置10、寝台装置20およびコンソール装置30の動作を制御することによって、X線検査装置1の全体の制御を行う処理部である。具体的には、システム制御部36は、スキャン制御部33を制御することにより、架台装置10および寝台装置20による被検体40のスペクトルのデータの収集動作を制御する。また、システム制御部36は、画像処理部34を制御することにより、断面画像の再構成処理を制御する。また、システム制御部36は、画像記憶部35から断面画像を読み出して、表示装置32に断面画像を表示させる。
なお、データ収集部16によって、収集したスペクトルのデータから所定のエネルギー帯ごとのサイノグラムが生成されるものとしているが、これに限定されるものではない。すなわち、データ収集部16は、収集したスペクトルのデータを画像処理部34に送信し、画像処理部34によって、スペクトルのデータから所定幅のエネルギー帯ごとのサイノグラムが生成されるものとしてもよい。
図2は、サイノグラムを説明する図である。図3は、検出器の特定のチャネルで検出されたエネルギーのスペクトルの一例を示す図である。図4は、被検体サイノグラムの例を示す図である。図5は、検出器に入射するX線の物理現象を説明する図である。図6は、検出スペクトルを説明する図である。図2〜6を参照しながら、サイノグラム、X線のエネルギーのスペクトル、および検出器13で検出される検出スペクトル等について説明する。
架台装置10のデータ収集部16は、検出器13により検出された、図3に示すようなエネルギーごとのフォトン数で示されるスペクトル(検出スペクトル)からサイノグラムを生成する。ここで、サイノグラムとは、図2に示すサイノグラム1001のように、X線管11のビュー毎、検出器13のチャネル毎の測定値を並べたデータであり、以降にお説明では、ビュー毎かつチャネル毎の測定値を画素値と見立て、サイノグラムを画像として扱うものとする。このうち、X線管11から照射されたX線が被検体40を透過して検出器13により検出されたスペクトル(図3参照)から生成されたサイノグラムを被検体サイノグラムというものとする。そして、被検体40を配置せずに、X線が空気だけを通過して検出器13により検出されたスペクトルから生成されたサイノグラムを空気サイノグラムというものとする。被検体サイノグラムおよび空気サイノグラムの画素値は、例えば、検出器13により測定値として検出されたフォトン数である。
また、検出器13は、ビュー毎、かつ、チャネル毎にエネルギーごとのフォトン数で示されるスペクトルを検出するので、データ収集部16は、X線管11の1周分のX線スキャンにより、図4に示すような、エネルギーごとの被検体サイノグラム1011を得ることができる。図4に示す例では、スペクトルを4つのエネルギー帯に分割して、エネルギー帯ごとに4つの被検体サイノグラム1011a〜1011dが得られる場合を示している。なお、図4では、4つのエネルギー帯に分割する例を示したが、この分割数に限定されるものではない。また、復元画像(減弱係数を示す画像)および後述する密度画像のS/N比の向上の観点から、再構成時、および物質密度の推定時に使用するエネルギー帯は、フォトン数が均一であることが望ましい場合がある。これを実現するために、例えば、以下の2の方法がある。
(方法1)後述の第1投影データを作成する段階で、フォトン数が均一となるようなエネルギー帯で分割する。
(方法2)まず、細かく分割しておき(例えば、1[keV]毎に分割)、再構成または物質密度の推定を行う段階で、フォトン数を合算する。
なお、被検体サイノグラムは、画素によってスペクトル形状が異なるので、空気サイノグラムに基づいて分割するものとしてもよい。
次に、図5を参照しながら、検出器13に入射するX線に対して発生する物理現象について説明する。検出器13は、図5に示すように、シンチレータ50と、接着層51と、SiPM(Si Photomultiplier)52と、を有するいわゆる間接変換型の検出器であるものとする。
シンチレータ50は、X線ビーム11b(図1参照)が入射すると、X線よりも波長の長い紫外線、可視光線または赤外線のうち少なくともいずれかを含む電磁波(以下、シンチレーション光というものとする)に変換する部材である。シンチレータ50は、反射板50aによって、回転フレーム12の周方向、および、スライス方向にマットリックス状に画定する板部材である。反射板50aは、X線は透過させるが、シンチレーション光を反射する性質を有する。
接着層51は、シンチレータ50の出射面側と、SiPM52の入射面側とを接着する層である。
SiPM52は、アバランシェフォトダイオード(Avalanche Photodiode:APD)とクエンチ抵抗との直列回路を複数並列に接続して構成された光電変換素子である。SiPM52は、シンチレータ50によりX線からシンチレーション光が接着層51を透過して入射すると、シンチレーション光から電流に光電変換する。シンチレータ50が反射板50aにより画定された各セルに対向するSiPM52の部分が、検出器13の各チャネルを形成する。
検出器13における特定のチャネルにX線が入射すると、シンチレータ50を構成する元素と反応する光電変換が生じ、X線のエネルギーのうち一部が蛍光X線として再放射され、残りのエネルギーはシンチレータ50と反応してシンチレーション光(可視光)となる。そして、蛍光が、図5に示すように、反射板50aを透過して他のチャネル(隣接するチャネル等)に移動する場合がある。この現象を、後述の蛍光と区別するために「エスケープ」というものとする。エスケープが発生すると、他のチャネルに移動したX線は、当該特定のチャネルにおけるフォトンの検出に寄与しないことになる。
また、検出器13における特定のチャネルの周囲のチャネルにX線が入射すると、上述と同様に、シンチレータ50を構成する元素と反応して光電効果を生じ蛍光を発する。そして、蛍光が、図5に示すように、周囲のチャネルから反射板50aを透過して、当該特定のチャネルに移動する場合がある。この現象を、上述のエスケープと区別して、単に「蛍光」というものとする。蛍光が発生すると、他のチャネルに入射したX線のエネルギーの一部が、当該特定のチャネルで検出されることになる。
また、シンチレータ50に入射したX線は、上述の光電変換だけではなく、シンチレータ50を構成する元素に衝突することによるコンプトン散乱およびレイリー散乱(以下、単に「散乱」という)が生じて、これにより進行方向が屈曲したX線が、図5に示すように、周囲のチャネルから反射板50aを透過して、当該特定のチャネルに移動する場合がある。この散乱が生じると、他のチャネルに入射したX線のエネルギーの一部またはすべてが、当該特定のチャネルで検出されることになる。
さらに、シンチレータ50に入射したX線のエネルギーの一部またはすべてが、シンチレータ50と反応してシンチレーション光に変換される。シンチレーション光は、上述のように、反射板50aを反射するが、図5に示すように、反射板50aが接着層51にまで延出して構成されていない場合、当該特定のチャネルの周囲のチャネルで発生したシンチレーション光は、接着層51を透過して、当該特定のチャネルのSiPM52で検出される場合がある。この現象を、以下では「クロストーク」というものとする。クロストークが発生すると、本来、他のチャネルに入射したX線であって、そのチャネルでフォトンが検出されるべきはずのX線のエネルギーの一部が、当該特定のチャネルで検出されることになる。
以上のようなエスケープ、蛍光、散乱およびクロストークの発生の影響によって、図6に示すように、被検体40を透過して検出器13の特定のチャネルに入射したX線のスペクトルである入射スペクトル2001は、検出器13によって、検出スペクトル2011のような波形のスペクトルで検出される。検出スペクトル2011は、図6に示すように、当該特定のチャネルに入射したX線のうちすべてのエネルギーが当該特定のチャネルに対応するSiPM52で検出された分のスペクトルと、スペクトル2012〜2014との和で表される。スペクトル2012は、エスケープにより、エネルギーの一部が他のチャネルに移動し、残りのエネルギーが当該特定のチャネルで検出されたスペクトルである。スペクトル2013は、蛍光により、他のチャネルに入射したX線のエネルギーの一部が当該特定のチャネルで検出されたスペクトルである。スペクトル2014は、クロストークおよび散乱により、他のチャネルに入射したX線のエネルギーの一部またはすべてが当該特定のチャネルで検出されたスペクトルである。すなわち、検出スペクトル2011は、本来、検出器13において検出されるべき正しい波形である入射スペクトル2001から歪んだ波形のスペクトルとなる。本実施形態に係るX線検査装置1の画像処理部34は、この歪んだ波形である検出スペクトル2011を、後述の補正動作により、入射スペクトル2001の波形に近づける動作を実行する。以下に、この画像処理部34の構成および動作について詳述する。
なお、検出器13は、図5に示すように、いわゆる間接変換型の検出器であるものとして説明したが、これに限定されるものではなく、シンチレータがなく、直接、X線を検出する直接変換型の検出器であってもよい。直接変換型の検出器であっても、検出器を構成する素材内をX線が透過する際に、上述のエスケープ、蛍光、クロストークおよび散乱と同等の作用を有する物理現象が起こり得る。
図7は、第1の実施形態の画像処理部のブロック構成の一例を示す図である。図7を参照しながら、本実施形態の画像処理部34のブロック構成について説明する。
図7に示すように、画像処理部34は、記憶部341と、生成部342と、再構成部343と、を備えている。
記憶部341は、後述する検出器応答データ(応答情報)を記憶する機能部である。なお、記憶部341は、画像処理部34が備える構成に限定されるものではなく、例えば、図1に示す画像記憶部35が記憶部341の機能を担うものとしてもよい。
生成部342は、データ収集部16から、被検体40のサイノグラムである被検体サイノグラムを第1投影データとして受信し、記憶部341から検出器応答データを読み出し、第1投影データおよび検出器応答データに基づいて、サイノグラムの形式である第2投影データを生成する機能部である。ここで、第1投影データは、検出器13によりビュー毎、かつ、チャネル毎に検出された、エネルギーごとのフォトン数で示されるスペクトルに基づいて、データ収集部16により生成された所定幅のエネルギー帯ごとの被検体サイノグラム(第1サイノグラム)の集合体で構成されるデータである。第2投影データも、第1投影データと同様のデータ形式であり、所定幅のエネルギー帯ごとの被検体サイノグラム(第2サイノグラム)の集合体で構成されるデータである。以下、第1投影データおよび第2投影データと同様のデータ形式であるデータを、単に「投影データ」という場合があるものとする。以下の説明では、第1投影データおよび第2投影データは、0〜140[keV]までのエネルギーにおけるフォトン数のスペクトルに関して、1[keV]毎のフォトン数、すなわち、0〜1[keV]、1〜2[keV]、・・・、139〜140[keV]のそれぞれのエネルギー帯(第1エネルギー帯、第2エネルギー帯)におけるフォトン数を画素値とする140個の被検体サイノグラムの集合体であるものとして説明する。
生成部342は、具体的には、下記の式(1)によって、第1投影データおよび検出器応答データから、第1投影データに含まれる歪みを補正した第2投影データを算出して生成する。
式(1)におけるpは、第1投影データおよび第2投影データである被検体サイノグラムを構成する画素を示す。ypは、生成部342による演算処理の対象となる第1投影データを示すデータであり、具体的には、第1投影データを構成するそれぞれの被検体サイノグラムの画素pの画素値(フォトン数)で構成される140次元のベクトルである。xpは、生成部342によって生成される第2投影データを示すデータであり、具体的には、第2投影データを構成するそれぞれの被検体サイノグラムの画素pの画素値(フォトン数)で構成される140次元のベクトルである。Mは、検出器13の特定のチャネル(画素pに対応するチャネル)に入射したX線が当該特定のチャネルで検出されるエネルギー応答を示し、すべてのエネルギーが当該特定のチャネルで検出された分のスペクトル、エスケープにより当該特定チャネルで検出された分のスペクトル、および、散乱により当該特定のチャネルで検出された分のスペクトル等から成り、140×140の行列である。Mのj列の成分は、検出器13の特定のチャネルに入射したX線のエネルギーのうち、(j−1)〜j[keV]のエネルギー帯のX線に対する応答スペクトルを示す。
また、式(1)におけるxiは、第2投影データを構成するそれぞれの被検体サイノグラムの画素pの周囲の画素iの画素値(フォトン数)で構成される140次元のベクトルである。Ciは、検出器13の特定のチャネル(画素pに対応するチャネル)の周囲の画素iに対応するチャネルに入力したX線のうち、蛍光、クロストークおよび散乱により当対特定のチャネルで検出された分のスペクトルに対する応答を示すデータであり、140×140の行列である。Ciのj列の成分は、検出器13の画素iに対応するチャネルに入射したX線のエネルギーのうち、(j−1)〜j[keV]のエネルギー帯であって、特定のチャネル(画素pに対応するチャネル)でX線が検出された場合の応答スペクトルを示す。Nは、画素pと同一のビューにおける、画素pの周囲の所定範囲(チャネル方向およびスライス方向の範囲)の画素の集合を示し、上述の画素iは集合Nに含まれる。
以下の説明においては、行列Mおよび行列Ciは、第1投影データおよび第2投影データである被検体サイノグラムのすべて画素で不変であるものとして説明する。ただし、行列Mおよび行列Ciは、検出器13の製造誤差に起因する特定のばらつき等に応じて、被検体サイノグラムの画素毎に設定するものとしてもよい。例えば、検出器13の特定のばらつきをチャネル毎に応答特性を計測しておき、行列Mおよび行列Ciを設定するものとすればよい。
ここで、検出器応答データとは、X線が入射した検出器13の応答特性として出力するスペクトルに含まれる誤差に寄与する物理現象の程度の大きさを示すデータである。検出器応答データは、一例として、検出器13の特定のチャネルに入射するX線のエネルギー毎の、当該特定のチャネルで生じるエスケープの発生確率、蛍光、クロストークおよび散乱等により周囲のチャネルから当該特定のチャネルで検出されるX線のエネルギーのスペクトルの情報、ならびに、検出されるエネルギーのばらつき等のうち少なくともいずれかに基づくデータである。具体的には、検出器応答データは、上述の式(1)における行列Mおよび行列Ciのデータを示す。すなわち、検出器応答データとは、上述のように、第1投影データを補正して第2投影データを生成するためのデータを示すものとする。
生成部342は、上述の式(1)に示すように、検出器応答データによるスペクトルの変化を考慮した結果が、第1投影データに近くなる投影データを第2投影データとして生成していることになる。
なお、上述の式(1)では、検出器応答データによるスペクトルの変化を考慮した結果と第1投影データの近さを、ベクトルのL2ノルムで定義したが、近さの算出方法はこれに限定されず、ベクトル間で定義される任意の距離を用いることができる。例えば、L1ノルム、L2ノルム、またはL∞ノルム等のLpノルムを用いることができる。また、ノルムの算出時は、成分ごとに重みを付与してもよい。また、ベクトルの成分間の相関を考慮した距離を用いてもよい。例えば、マハラノビス距離を用いることができる。
また、第1投影データおよび第2投影データにおけるエネルギー帯の幅は、上述した1[keV]に限定されるものではなく、エネルギーの単位も[keV]に限定されるものではなく、また、スペクトルのエネルギーの範囲も0〜140[keV]に限定されるものではない。すなわち、異なる幅のエネルギー帯であってもよく、エネルギーの単位が[J]または[cal]等の異なる単位であってもよい。また、第1投影データおよび第2投影データにおける単位はエネルギーの単位に限定されるものではない。例えば、第1投影データにおける単位は、検出器13から読み出したデジタル信号値を単位とし、第2投影データにおける単位は、[keV]とするものとしてもよい。
また、第1投影データおよび第2投影データは、0〜140[keV]までのエネルギーにおけるフォトン数のスペクトルに関して、140個のエネルギー帯におけるそれぞれのフォトン数を画素値とする140個の被検体サイノグラムの集合体としたが、第1投影データと第2投影データとで、エネルギー帯の数が同一であることに限定されない。例えば、第1投影データについては、0〜140[keV]までのエネルギーについて、5[keV]毎のエネルギー帯(第1エネルギー帯)のフォトン数を画素値とする28個の被検体サイノグラムの集合体とし、第2投影データについては、上述のように、1[keV]毎のエネルギー帯(第2エネルギー帯)のフォトン数を画素値とする140個の被検体サイノグラムの集合体としてもよい。この場合、第1投影データについてのエネルギー帯の数と、第2投影データについてのエネルギー帯の数との差は、行列Mおよび行列Ciの次元を28×140とすればよい。この時、行列Mおよび行列Ciのj行目の成分は、各行列を140×140の次元で作成したときの、(5×j+1)行目から(5×j+5)行目の成分を、各列ごとに加算した値とすればよい。
また、生成部342は、実際に上述の式(1)によって第2投影データを算出する場合、下記の式(2)に示すような反復計算により算出する。
生成部342は、まず、式(2)において、第2投影データの初期値を任意に設定する。例えば、生成部342は、第2投影データの初期値を第1投影データとするものとしてもよい。そして、生成部342は、xp (k+1)を(k+1)回目の反復計算後の第2投影データとし、このxp (k+1)を次の更新のための式(2)の計算にxi (k)として利用する。なお、反復計算の回数は、例えば、所定回数行うものとしてもよく、式(2)により計算した第2投影データが、その直前に式(2)により計算した第2投影データとの差が所定の閾値以下となった場合に反復計算を停止するものとしてもよい。
なお、上述の検出器応答データは、単位時間あたりに検出器13に入射するX線のフォトンの数によって変化する場合がある。例えば、X線のフォトンが連続して検出器13に入射した場合、検出器13の反応が不十分になり、見かけ上、実際よりも低いエネルギーのX線のフォトン数として検出される。そこで、生成部342は、検出器応答データを、第1投影データのフォトン数に応じて切り替えて、第2投影データを生成するものとしてもよい。この場合、例えば、記憶部341には、単位時間あたりに検出器13に入射するフォトン数を10段階に分け、それぞれに対する検出器応答データを記憶しておく。そして、生成部342は、第1投影データの被検体サイノグラムの画素毎にそのフォトン数を計算し、そのフォトン数に該当する検出器応答データを記憶部341から読み込んで用いる。また、検出器応答データは、第1投影データの被検体サイノグラムの画素のフォトン数に近い検出器応答データを合成して読み込むこともできる。例えば、フォトン数「500」に対する検出器応答データであるM500およびCi,500と、フォトン数「1000」に対する検出器応答データであるM1000およびCi,1000が記憶されているとする。この時、第1投影データの被検体サイノグラムの画素のフォトン数が「700」であった場合は、フォトン数「500」および「1000」それぞれへのフォトン数「700」の近さに応じて検出器応答データを加重平均する。すなわち、検出器応答データとして、(3/5)M500+(2/5)M1000および、(3/5)Ci,500+(2/5)Ci,1000を計算して読み込む。また、記憶部341に、フォトン数の変化による検出器応答データの変動分を補正するための係数を記憶しておき、第1投影データの被検体サイノグラムの画素のフォトン数に応じて係数を調整し、検出器応答データを補正した後に読み込むこともできる。なお、被検体40を透過し検出器13の特定のチャネルに入力するX線の入射スペクトルが、当該特定のチャネルの空間的および時間的近傍のチャネルに入射するX線の入射スペクトルと十分に近似しているとして、近傍のデータを合成することで第1投影データの値とすることもできる。合成の一例としては、加重平均を用いることができる。ここで、第1投影データのサイノグラムの当該画素に対して、空間的近傍のチャネルとは、チャネル方向及びスライス方向に近傍のチャネルであり、時間的近傍のチャネルとは、近傍のビューのデータ、あるいは、検出器を2回転以上回転させて収集した投影データにおける、ちょうど1回転分ビューが異なるビューの近傍のビューのデータである。これにより、線量が低く、第1投影データがフォトンゆらぎの影響により、精度が高く観測できない場合にも、高精度に第2投影データを生成可能である。
再構成部343は、生成部342により生成された第2投影データに含まれる被検体サイノグラムのうち、復元の対象とするエネルギー帯の被検体サイノグラムを再構成して復元画像を生成する機能部である。復元画像は、例えば、線減弱係数を画素値とする画像である。再構成の方法としては、例えば、逆投影法の一例であるフィルタ補正逆投影法(FBP(Filtered Back Projection)法)を用いることができる。ここでは、第2投影データの被検体サイノグラムを再構成する方法として、FBP法を採用した場合について説明する。この場合、再構成部343は、被検体40を配置せずに、X線が空気だけを通過して検出器13により検出されたスペクトルにおいて上述の復元の対象とするエネルギー帯に対応する部分から生成された空気サイノグラムである基準データI0を有しているものとする。また、第2投影データに含まれる被検体サイノグラムのうち、復元の対象となるエネルギー帯の被検体サイノグラムを復元対象データIとする。復元対象データIは、復元の対象となるエネルギー帯に属するエネルギーに関して、第2投影データのフォトン数の和により算出する。
まず、再構成部343は、まず、復元対象データIおよび基準データI
0から、下記の式(3)により線減弱係数の積分値Mを算出する。
式(3)の(m,n)はm番目のチャネル、かつ、n番目のビューのデータであることを示す。式(3)に示す積分値M(m,n)は、n番目のビューにおいて、X線管11から照射されたX線が検出器13のm番目のチャネルに到達するまでに通過する経路に沿って被検体40の線減弱係数を積分した値となる。
次に、再構成部343は、算出した積分値M(m,n)に対して、チャネル方向に1次元フーリエ変換をする。次に、再構成部343は、1次元フーリエ変換をした値に対して、ランプフィルタまたはShepp−Loganフィルタ等の高域通過フィルタによって周波数方向にフィルタ処理をして、1次元逆フーリエ変換をする。そして、再構成部343は、1次元逆フーリエ変換をしたデータを、すべてのビュー毎に逆投影して加算することによって、再構成された再構成画像である復元画像を生成する。
なお、再構成部343は、逆投影法の一例であるFBP法を用いて再構成する場合を示したが、再構成する方法はこれに限定されるものではなく、逐次近似法等の各種の再構成方式を用いてもよい。例えば、逐次近似法では、疑似的に仮の画像を予め用意し、各ビューにおいてX線を被検体に照射して減弱率を算出する方法である。そして、仮の画像で算出した減弱率が、実際に検出器13において検出された測定値(減弱率)より小さい場合、仮の画像の画素値を増加させていく。逆に、仮の画像で算出した減弱率が、実際に検出器13において検出された測定値より大きい場合、仮の画像の画素値を減少させていく。この動作を繰り返すことによって、仮の画像で算出される減弱率が、実際に検出器13で検出された測定値(減弱率)と等しくなるように変更して再構成画像を得る。逐次近似法には、ART(Algebraic Reconstruction Technique)法、OS−EM(Ordered Subset Expectation Maximization)法、およびML−EM(Maximum Likelihood Expectation Maximization)法等の種々の方法がある。
図8は、第1の実施形態の画像処理部の動作の一例を示すフローチャートである。図8を参照しながら、本実施形態の画像処理部34による画像処理動作について説明する。
<ステップS11>
画像処理部34の生成部342は、データ収集部16から、被検体40のサイノグラムである被検体サイノグラムを第1投影データとして受信して取得する。そして、ステップS12へ移行する。
<ステップS12>
生成部342は、画像処理部34の記憶部341から検出器応答データを読み出して取得する。そして、ステップS13へ移行する。
<ステップS13>
生成部342は、上述の式(1)を用いて、取得した第1投影データおよび検出器応答データから、第1投影データに含まれる歪みを補正した第2投影データを生成する。この場合、検出器応答データは、上述の式(1)における行列Mおよび行列Ciのデータを示す。生成部342は、上述の式(1)に示すように、検出器応答データを作用させた投影データが、第1投影データに近くなる投影データを第2投影データとして生成することになる。生成部342は、生成した第2投影データを、画像処理部34の再構成部343へ送る。そして、ステップS14へ移行する。
<ステップS14>
再構成部343は、生成部342により生成された第2投影データに含まれる被検体サイノグラムのうち、復元の対象とするエネルギー帯の被検体サイノグラムを再構成して復元画像を生成する。具体的には、再構成部343は、まず、第2投影データに含まれる被検体サイノグラムのうち、復元の対象となるエネルギー帯の被検体サイノグラムである復元対象データI、および基準データI0から、上述の式(3)により線減弱係数の積分値M(m,n)を算出する。次に、再構成部343は、算出した積分値M(m,n)に対して、チャネル方向に1次元フーリエ変換をする。次に、再構成部343は、1次元フーリエ変換をした値に対して、ランプフィルタまたはShepp−Loganフィルタ等の高域通過フィルタによって周波数方向にフィルタ処理をして、1次元逆フーリエ変換をする。そして、再構成部343は、1次元逆フーリエ変換をしたデータを、すべてのビュー毎に逆投影して加算することによって、再構成された再構成画像である復元画像を生成する。
以上のステップS11〜S14の動作によって、画像処理部34による画像処理が実行される。
以上のように、検出器13の特定のチャネルに入射するX線のエネルギー毎の、当該特定のチャネルで生じるエスケープの確率、ならびに、蛍光、クロストークおよび散乱等により周囲のチャネルから当該特定のチャネルで検出されるX線のエネルギーのスペクトルの情報等である検出器応答データを用いて、生成部342は、第1投影データに含まれる歪みを補正した第2投影データを算出するものとしている。すなわち、生成部342は、上述の式(1)に示すように、検出器応答データを作用させた投影データが、第1投影データに近くなる投影データを第2投影データとして生成するものとしている。そして、再構成部343は、生成部342により生成された第2投影データに含まれる被検体サイノグラムのうち、復元の対象とするエネルギー帯の被検体サイノグラムを再構成して復元画像を生成する。これによって、被検体の組成によらず、所定のエネルギー帯ごとの被検体サイノグラムの画素値を理論値(第1投影データ)に近づけるように補正でき、検出器13により検出されたX線のスペクトルの歪みの補正の精度を向上させ、最終的に再構成された復元画像の正確性を向上させることができる。例えば、第2投影データの各被検体サイノグラムからスペクトルを復元したものとした場合の復元スペクトル(以下、単に「第2投影データの復元スペクトル」という)は、上述の被検体40を透過し検出器13に入射したX線のスペクトル(入射スペクトル)の形状に近づけることができる。
なお、検出器13は、回転フレーム12の周方向に並んだチャネル(検出素子)ごとにエネルギーごとのフォトン数で示されるスペクトルを検出するものとしたが、上述のように、検出器13は、被検体40の体軸方向にも検出素子が配列されている。したがって、体軸方向(スライス方向)における検出素子のリング状の配列ごとにサイノグラムが生成され、上述の画像処理が実行されるものとしてよい。または、ヘリカルスキャンのように天板22を移動させながら回転フレーム12を連続回転させる場合においては、同一の周方向のチャネル(検出素子)により検出されたデータのみ使用するのではなく、体軸方向(スライス方向)にずれたチャネルによって検出されたデータにより補間してサイノグラムが生成されるものとしてもよい。また、デュアルエナジーX線CT装置のように、X線管11から照射されるX線のエネルギーを2種類に分けて1周ごとに切り替えて照射(例えば、1周目は140[keV]、2周目は80[keV])して、異なるエネルギーのスペクトルを合成して、サイノグラムを生成するものとしてもよい。
また、再構成部343によって生成される復元画像の画素値は線減弱係数であるものとしたが、これに限定されるものではなく、CT値等のX線の減弱の多寡を表す値であれば、いずれも有効である。同様に、サイノグラムの画素値も、X線の量を示す値であればよい。例えば、X線の量そのもの、もしくはフォトン数の多寡を表す値、または、X線の量もしくはフォトン数の変化率を表す値であってもよい。
<変形例1>
本変形例の画像処理部34aについて、第1の実施形態の画像処理部34と相違する点を中心に説明する。本変形例の画像処理部34aは、図7に示す第1の実施形態の画像処理部34が備える記憶部341および生成部342の代わりに、それぞれ、後述する図9に示す記憶部341aおよび生成部342aを備えている。
図9は、第1の実施形態の変形例1の画像処理部のブロック構成の一例を示す図である。図10は、検出スペクトルが補正された復元スペクトルの例を示す図である。図9を参照しながら、本変形例の画像処理部34aのブロック構成および動作について説明する。
図9に示すように、画像処理部34aは、記憶部341aと、生成部342aと、再構成部343と、を備えている。
記憶部341aは、後述するように、第1投影データを第2投影データに変換するためのフィルタ係数を記憶する機能部である。
生成部342aは、データ収集部16から、被検体40のサイノグラムである被検体サイノグラムを第1投影データとして受信し、記憶部341aから、第1投影データを第2投影データに変換するためのフィルタ係数を読み出し、第1投影データおよびフィルタ係数に基づいて、サイノグラムの形式である第2投影データを生成する機能部である。ここで、例えば、被検体40を透過し検出器13の特定のチャネルに入力するX線の入射スペクトルは、当該特定のチャネルの空間的近傍のチャネルに入射するX線の入射スペクトルと十分に近似しているものとする。この場合、上述の式(1)は、下記の式(4)に置き換えることができる。
ただし、Hは、下記の式(5)で表される140×140の行列である。
生成部342aは、上述の式(4)によって、第1投影データ、および行列Hのデータ(上述の、検出器応答データに相当)から、第1投影データに含まれる歪みを補正した第2投影データを算出することも可能である。ただし、ここでは、ノイズまたは検出器13の検出誤差への耐性を持たせるために、下記の式(6)のように、式(4)に示す値に正則化項λf(連続性を示す情報)を加算するものとする。
式(6)におけるλは、正則化の強さを表す重みであり、f(x
p)は、x
pのエネルギー方向のスペクトルの滑らかさ(連続性)を表す関数である。関数f(x)は、例えば、下記の式(7)で表される。
あるいは、関数f(x)は、下記の式(8)のように表されるものとしてもよい。
式(7)および(8)におけるx(j)は、xのj番目の成分(フォトン数)である。式(7)および(8)に示すf(x)は、いずれも、隣接する成分(フォトン数)の差が大きいほど、値が大きくなる関数である。
上述の式(6)は、勾配法または劣勾配法等の反復アルゴリズムで解くことが可能である。また、関数f(x)が上述の式(7)で示される場合、解析的に解くことも可能であり、以下では、この場合について説明する。
関数f(x)が上述の式(7)で表される場合、上述の式(6)は、下記の式(9)と等価である。
式(9)におけるSは、下記の式(10)で表される(140+139)×140の行列である。
また、式(9)におけるy’
pは、下記の式(11)で表される(140+139)次元のベクトルである。
上述の式(9)の解であるx
pは、行列Sの擬似逆行列S
+を用いて、下記の式(12)によって算出できる。また、擬似逆行列S
+は、行列Sの特異値分解を用いて算出可能である。
ここで、上述した記憶部341aが記憶するフィルタ係数とは、擬似逆行列S+の成分の値である。生成部342aは、この式(12)によって、第1投影データに対して、記憶部341aから読み出したフィルタ係数を畳み込むことによって、第1投影データに含まれる歪みを補正した第2投影データを生成する。すなわち、記憶部341aが記憶するフィルタ係数は、第1投影データを補正して第2投影データを生成するためのデータであるので、上述の検出器応答データと捉えることができる。
再構成部343は、生成部342aにより生成された第2投影データに含まれる被検体サイノグラムのうち、復元の対象とするエネルギー帯の被検体サイノグラムを再構成して復元画像を生成する機能部である。再構成部343による再構成の方法は、第1の実施形態と同様である。
ここで、被検体40を透過して検出器13の特定のチャネルに、図10に示す入射スペクトル2101で示されるX線が入射されたものとする。この場合、検出器13によって検出されるスペクトルは、上述したように、エスケープ、蛍光、クロストークおよび散乱等によって、入射スペクトル2101と比較して歪んだ検出スペクトル2111となる。そして、データ収集部16は、検出器13から収集したスペクトル(検出スペクトル2111はその一例)のデータそれぞれに対して増幅処理またはA/D変換処理等を行なって、エネルギー帯ごとの被検体サイノグラムを生成し、第1投影データとして生成部342aに送る。さらに、生成部342aによって、上述の式(6)の正則化項が加味された式(12)を用いて、算出された第2投影データの各被検体サイノグラムからスペクトルを復元した場合のスペクトルが、図10に示す復元スペクトル2121である。このように、生成部342aによって、正則化項が加味された式(12)を用いて第2投影データが算出されることにより、復元スペクトル2121は、検出スペクトル2111と比較して、入射スペクトル2101に近似した形状に近づけることができる。また、生成部342aは、単純化された式(12)を用いて、第2投影データを生成することができるので、演算負荷を低減することができる。
なお、記憶部341aは、上述の式(12)に示す擬似行列S+の成分の値を係数として記憶するものとし、生成部342aは、この式(12)によって、第1投影データに対して、記憶部341aから読み出したフィルタ係数を畳み込むことによって、第2投影データを生成するものとしたが、これに限定されるものではない。例えば、記憶部341aは、擬似行列S+の成分の値である係数を記憶するのではなく、検出器応答データとしての行列Hおよび重みλを記憶しておくものとしてもよい。この場合、生成部342aは、記憶部341aに記憶された検出器応答データを用いて、式(6)〜(12)により、第2投影データを算出して生成するものとしてもよい。
<変形例2>
本変形例の再構成部343の動作について、第1の実施形態の再構成部343と相違する点を中心に説明する。本変形例の再構成部343は、再構成した復元画像の画素である線減弱係数、および特定の物質の質量減弱係数に基づいて、その特定の物質の密度を推定する。以下、その物質密度の推定動作の詳細を説明する。
本変形例の画像処理部34の再構成部343は、上述の第1の実施形態のように、第2投影データから、対象とする1つのエネルギー帯の被検体サイノグラムを再構成するのではなく、複数のエネルギー帯の被検体サイノグラムを再構成して、画素値が線減弱係数である複数の復元画像を生成する。ここで、例えば、被検体40に含まれる物質として、水およびヨードを想定し、それぞれの物質密度を推定する場合について説明する。
まず、再構成部343は、第2投影データに含まれる被検体サイノグラムのうち、特定のエネルギー帯(例えば、35〜50[keV])および特定のエネルギー帯とは異なるエネルギー帯(例えば、55〜70[keV])のそれぞれ(第2エネルギー帯)の被検体サイノグラムを再構成して、画素値が線減弱係数である2つの復元画像を生成する。そして、再構成部343は、下記の式(13)に示す連立方程式から、復元画像の画素毎に物質密度を算出する。
式(13)の(s,t)は、復元画像の座標を示す。ρwおよびρIは、それぞれ水およびヨードの物質密度であり、式(13)の算出対象となる値である。vw,1およびvI,1は、それぞれ特定のエネルギー帯の水の質量減弱係数およびヨードの質量減弱係数であり、既知の値である。vw,2およびvI,2は、それぞれ特定のエネルギー帯とは異なるエネルギー帯の水の質量減弱係数およびヨードの質量減弱係数であり、既知の値である。μ1は、特定のエネルギー帯の被検体サイノグラムから再構成部343により再構成された復元画像の線減弱係数である画素値である。そして、μ2は、特定のエネルギー帯とは異なるエネルギー帯の被検体サイノグラムから再構成部343により再構成された復元画像の線減弱係数である画素値である。上述のように、vw,1、vI,1、vw,2およびvI,2は、既知の値であり、μ1およびμ2は、得られている値なので、式(13)は、ρwおよびρIについての連立方程式である。再構成部343は、ρwおよびρIについての連立方程式である式(13)を解くことにより、復元画像の画素(座標(s,t))毎に、ρwおよびρIを算出することができる。そして、再構成部343は、復元画像の各画素の画素値を、その画素に対応する物質密度ρwまたはρIに置換することによって、水またはヨードの密度画像を生成することができる。
以上のように、再構成部343は、正確性が向上した線減弱係数で構成される復元画像を生成するので、その復元画像から被検体40の特定の物質について正確な密度を示す密度画像を生成することができる。
なお、上述では、水およびヨードの密度を求める例を示したが、密度を求める物質はこれらに限定されるものではなく、例えば、骨またはガドリニウム等の任意の物質の密度画像を生成することができる。また、上述の例では、水およびヨードの2種類の物質の密度を求める例を示したが、これに限定されるものではなく、3種類以上の物質の密度を求めることも可能である。この場合、密度を計算する物質の種類の数と同数以上の種類のエネルギー帯の被検体サイノグラムから、画素値が線減弱係数である復元画像を生成するものとすればよい。
なお、上述した、再構成画像を生成した後に物質密度の推定を行う以外に、サイノグラムの各画素に対して物質密度の推定を行った後に再構成処理を行うことでも、物質密度を画素値として有する断面画像を生成することもできる。
例として、水およびヨードの物質密度画像を算出する方法を説明する。まず、第2投影データに含まれる被検体サイノグラムのうち、特定のエネルギー帯(例えば、35〜50[keV])(以下、エネルギー帯Aという)、および、特定のエネルギー帯とは異なるエネルギー帯(例えば、55〜70[keV])(以下、エネルギー帯Bという)に関して、それぞれのエネルギー帯に属するエネルギーの被検体サイノグラムの値の和を計算する(それぞれを、I1およびI2とする)。同様に、エネルギー帯Aとエネルギー帯Bに対応する、空気サイノグラムの値の和も計算する(それぞれを、I0,1およびI0,2とする)。
次に、下記の式(14)により、それぞれのエネルギー帯毎に線減弱係数の積分値M1およびM2を計算する。
そして、下記の式(15)に示す連立方程式から、サイノグラムの画素毎に物質の透過距離を計算する。
ここで、lwおよびlIは、それぞれサイノグラムの画素(m,n)に関連するX線の透過経路上にある水およびヨードの距離であり、算出対象の値である。また、zw,1およびzI,1は、それぞれエネルギー帯Aの水およびヨードの線減弱係数であり、zw,2およびzI,2は、それぞれエネルギー帯Bの水およびヨードの線減弱係数である。次に、lwおよびlIをFBP法等で再構成処理した後に、水およびヨードの密度をそれぞれ掛け合わせることで、水の密度を示す断面画像と、およびヨードの密度を示す断面画像を得ることができる。上述した例では、水およびヨードを物質の例にしたが、物質はこれらに限定されない。また、3つ以上の物質密度を計算することもできる。この場合は、エネルギー帯を物質の数よりも多く作成し、連立方程式の式の数を物質の数以上にすればよい。
(第2の実施形態)
本実施形態の画像処理部34bについて、第1の実施形態の変形例1の画像処理部34aと相違する点を中心に説明する。本実施形態の画像処理部34bは、上述の式(6)における正則化項λfにおける、式(7)等で示される関数f(x)について、隣り合うエネルギー帯の画素値(フォトン数)の差に対して重みを加味して算出する。以下、重みを加味した算出動作の詳細を説明する。
図11は、第2の実施形態の画像処理部のブロック構成の一例を示す図である。図12は、X線管から出射された出射スペクトルおよび被検体スペクトルの一例を示す図である。図11を参照しながら、本実施形態の画像処理部34bのブロック構成について説明する。
図11に示すように、画像処理部34bは、記憶部341bと、生成部342bと、再構成部343と、算出部344と、を備えている。
記憶部341bは、上述の式(6)に示す行列Hおよび重みλを検出器応答データとして記憶する機能部である。
生成部342bは、データ収集部16から、被検体40のサイノグラムである被検体サイノグラムを第1投影データとして受信し、記憶部341bから検出器応答データを読み出し、第1投影データおよび検出器応答データに基づいて、サイノグラムの形式である第2投影データを生成する機能部である。また、第1の実施形態の変形例1と同様に、被検体40を透過し検出器13の特定のチャネルに入力するX線の入射スペクトルは、当該特定のチャネルの近傍のチャネルに入射するX線の入射スペクトルと十分に近似しているものとする。
上述の式(6)は、上述の式(4)に正則化項λfを加算したものであるが、本実施形態では、この正則化項λfのうち関数f(x)を、下記の式(16)に示すように、隣り合うエネルギー帯の画素値(フォトン数)の差ごとに第1重みw
jを乗算して算出する。
あるいは、関数f(x)は、下記の式(17)のように表されるものとしてもよい。
また、関数f(x)が上述の式(16)で表される場合、上述の式(6)は、下記の式(18)と等価である。
ただし、式(18)におけるWは、下記の式(19)で表される(140+139)×(140+139)の行列である。
上述の式(18)の解であるx
pは、(140+139)×140の行列である行列WSの擬似逆行列(WS)
+を用いて、下記の式(20)によって算出できる。また、擬似逆行列(WS)
+は、行列Sの特異値分解を用いて算出可能である。
生成部342bは、後述する算出部344により算出された第1重みwj、および、記憶部341bから読み出した検出器応答データ(行列H、重みλ)から、式(6)、(16)、および(18)〜(20)によって、第1投影データに含まれる歪みを補正した第2投影データを算出して生成する。
再構成部343は、生成部342bにより生成された第2投影データに含まれる被検体サイノグラムのうち、復元の対象とするエネルギー帯の被検体サイノグラムを再構成して復元画像を生成する機能部である。再構成部343による再構成の方法は、第1の実施形態と同様である。
算出部344は、生成部342bが第2投影データを算出するために、上述の式(16)の第1重みwjを算出する機能部である。算出部344は、例えば、X線管11から出射されるX線のスペクトルである出射スペクトルの情報、X線管11のターゲットの部材、X線管11のターゲットの角度、X線管11が備えるフィルタ、撮影時のX線管11の管電圧、使用した造影剤、および、想定される被検体40の組成のうちの少なくともいずれかの情報に基づいて、第2投影データの復元スペクトルがエネルギー方向に急峻に変換すると想定される程度を求めて第1重みwjを算出する。算出部344は、第1重みwjを、第2投影データの復元スペクトルにおいて、エネルギー方向に急峻に変化すると想定されるエネルギーでは小さくし、エネルギー方向に滑らかであると想定されるエネルギーでは大きくするように算出する。なお、算出部344は、上述のX線管11から出射されるX線のスペクトルである出射スペクトルの情報、X線管11のターゲットの部材、X線管11のターゲットの角度、X線管11が備えるフィルタ、撮影時のX線管11の管電圧、使用した造影剤、および、想定される被検体40の組成等の情報は、例えば、システム制御部36等から取得するものとすればよい。
以下、図12を参照しながら、算出部344による第1重みwjの算出の方法の具体例について説明する。算出部344は、例えば、図12(a)に示すように、X線管11から出射されるX線のスペクトルである出射スペクトル2200の情報を、システム制御部36等から取得する。算出部344は、出射スペクトル2200の形状から、エネルギーE1およびE2の部分でエネルギー方向に急峻に変化していると検出し、エネルギーE1およびE2の周辺のエネルギーで第1重みwjが小さくなるように算出する。算出部344は、例えば、エネルギーE1±3[keV]およびE2±3[keV]の範囲で、第1重みwjを「1」とし、それ以外のエネルギーでは第1重みwjを「10」とする。これは、出射スペクトル2200に含まれるエネルギー方向の急峻な変化は、被検体40の透過後にも、同じエネルギーにおいてエネルギー方向にスペクトルが急峻に変化することに基づいている。
なお、X線管11から出射されるX線の出射スペクトルは、その形状の情報を直接取得しなくても、X線管11のターゲットの部材、X線管11のターゲットの角度、X線管11が備えるフィルタ、および撮影時のX線管11の管電圧の情報等から出射スペクトルの形状の計算が可能であり、その計算結果を使用することもできる。また、出射スペクトルの形状を計算しなくても、X線管11のターゲットの部材、X線管11のターゲットの角度、X線管11が備えるフィルタ、および撮影時のX線管11の管電圧の情報等の少なくともいずれかの情報から、X線管11から出射されるX線の出射スペクトルが、エネルギー方向に急峻に変化するエネルギーの情報を得ることもできる。
また、被検体40にヨードまたはカドリニウム等の造影剤が含まれている場合、それぞれの造影剤のK吸収端のエネルギーで、被検体40の通過により減弱したX線のスペクトルである被検体スペクトル(検出器13に入射するという意味では入射スペクトルともいえる)がエネルギー方向で急峻に変化する。図12(b)は、造影剤としてヨードを含む被検体40をX線が減弱した場合のスペクトルの例として、被検体スペクトル2201を示している。図12(b)に示す被検体スペクトル2201は、エネルギーE3でスペクトルが急峻に変化している。したがって、算出部344は、被検体スペクトル2201の形状から、エネルギーE3の部分でエネルギー方向に急峻に変化していると検出し、エネルギーE3の周辺のエネルギーで第1重みwjが小さくなるように算出する。算出部344は、例えば、エネルギーE1±3[keV]およびE2±3[keV]の範囲に加えて、エネルギーE3±3[keV]における第1重みwjを「1」とし、それ以外のエネルギーでは第1重みwjを「10」とする。ここで、K吸収端のエネルギーは、造影剤ごとに固有のエネルギーであり、ヨードでは33[keV]、ガドリニウムでは50[keV]である。
また、造影剤の情報に基づく算出部344による第1重みwjの算出は、被検体40に造影剤が含まれるときだけ行われることが好ましいが、造影剤の有無が不明確な場合には、使用の可能性のある造影剤に対応させて第1重みwjの算出が行われるものとしてもよい。また、算出部344は、使用する造影剤の情報を、システム制御部36等から受けることもできる。また、ユーザが、入力装置31を介して、使用する造影剤の情報を入力してもよい。さらに、ユーザが、入力装置31を介して、第2投影データがエネルギー方向に急峻に変化すると想定されるエネルギーの値を、直接入力可能とするようにしてもよい。
以上のように、生成部342bは、算出部344により算出された第1重みwjを利用して、第2投影データを生成する。すなわち、生成部342bは、第1重みwjが大きいエネルギーではエネルギー方向の滑らかさを高めるように作用させ、第1重みwjが小さいエネルギーではエネルギー方向の滑らかさを高める作用を軽減するようにし、検出器応答データを作用させた投影データが、第1投影データに近くなる投影データを第2投影データとして生成する。
ここで、被検体40を透過して検出器13の特定のチャネルに、上述の図10に示す入射スペクトル2101で示されるX線が入射されたものとする。この場合、検出器13によって検出されるスペクトルは、上述したように、エスケープ、蛍光、クロストークおよび散乱等によって、入射スペクトル2101と比較して歪んだ検出スペクトル2111となる。そして、データ収集部16は、検出器13から収集したスペクトル(検出スペクトル2111はその一例)のデータそれぞれに対して増幅処理またはA/D変換処理等を行なって、エネルギー帯ごとの被検体サイノグラムを生成し、第1投影データとして生成部342bに送る。さらに、生成部342bによって、上述の式(6)の正則化項および第1重みwjが加味された式(20)を用いて、算出された第2投影データの各被検体サイノグラムからスペクトルを復元した場合のスペクトルが、図10に示す復元スペクトル2122である。このように、生成部342aによって、正則化項および第1重みwjが加味された式(20)を用いて第2投影データが算出されることにより、復元スペクトル2122は、検出スペクトル2111と比較して、入射スペクトル2101においてエネルギー方向に急峻に変化している部分についても近似した形状に近づけることができる。
なお、生成部342bは、さらに、スペクトルの誤差に対してエネルギー帯ごとに重み付けして、第2投影データを生成するものとしてもよい。このスペクトルの誤差に対してエネルギー帯ごとに重み付けする重みを、第2重みu
jとすると、上述の式(17)に示す行列Wの代わりに、下記の式(21)のように表される。
すなわち、生成部342bは、上述の第2重みuj、算出部344により算出された第1重みwj、および、記憶部341bから読み出した検出器応答データ(行列H、重みλ)から、式(6)、(16)、(17)、(19)および(21)によって、第2投影データを算出して生成するものとしてもよい。第2重みujは、第1投影データにおいて信頼性が低いエネルギー帯に対応するものを小さくすることが望ましい。例えば、低いエネルギーのスペクトルは、検出器13のノイズの影響を受けており、信頼性が低いといえる。そこで、例えば、所定のエネルギーより小さいエネルギー帯では、第2重みujを「1」とし、大きいエネルギー帯では、第2重みujを「10」とするものとすればよい。また、第1投影データにおいてフォトン数が少ないエネルギー帯についても、信頼性が低いと判断することができる。
なお、第2重みujは、被検体サイノグラムの画素毎に個別に算出してもよい。例えば、被検体40によるX線の減弱が少ない被検体サイノグラムの画素では、検出器13のノイズに対して影響が小さくなるため、小さいエネルギーにおいても画素値(フォトン数)の信頼度が高いため、第2重みujをすべてのエネルギーで均一にする。一方、被検体40によるX線の減弱が大きい被検体サイノグラムの画素では、検出器13のノイズの影響を考慮して、小さいエネルギーでは第2重みujを小さく設定するようにするといった、被検体サイノグラムの画素毎に第2重みujの調整を行うことができる。
図13は、第2の実施形態の画像処理部の動作の一例を示すフローチャートである。図13を参照しながら、本実施形態の画像処理部34bによる画像処理動作について説明する。
<ステップS21>
画像処理部34bの生成部342bは、データ収集部16から、被検体40のサイノグラムである被検体サイノグラムを第1投影データとして受信して取得する。そして、ステップS22へ移行する。
<ステップS22>
生成部342bは、画像処理部34bの記憶部341bから検出器応答データを読み出して取得する。そして、ステップS23へ移行する。
<ステップS23>
算出部344は、生成部342bが第2投影データを算出するために利用する上述の式(16)の第1重みwjを算出する。まず、算出部344は、例えば、X線管11から出射されるX線のスペクトルである出射スペクトルの情報、X線管11のターゲットの部材、X線管11のターゲットの角度、X線管11が備えるフィルタ、撮影時のX線管11の管電圧、使用した造影剤、および、想定される被検体40の組成のうちの少なくともいずれかの情報を、システム制御部36等から取得する。そして、算出部344は、取得した情報に基づいて、第2投影データの復元スペクトルがエネルギー方向に急峻に変換すると想定される程度を求めて第1重みwjを算出する。すなわち、算出部344は、第1重みwjを、第2投影データの復元スペクトルにおいて、エネルギー方向に急峻に変化すると想定されるエネルギーでは小さくし、エネルギー方向に滑らかであると想定されるエネルギーでは大きくするように算出する。算出部344は、算出した第1重みwjを、生成部342bに送る。そして、ステップS24へ移行する。
<ステップS24>
生成部342bは、算出部344から受け取った第1重みwj、および検出器応答データから、式(6)、(16)、および(18)〜(20)を用いて、第1投影データに含まれる歪みを補正した第2投影データを算出して生成する。この場合、検出器応答データは、上述の式(6)における行列Hおよび重みλのデータを示す。生成部342bは、検出器応答データを作用させた投影データが、第1投影データに近くなる投影データを第2投影データとして生成する。生成部342bは、生成した第2投影データを、画像処理部34の再構成部343へ送る。そして、ステップS25へ移行する。
<ステップS25>
再構成部343は、生成部342bにより生成された第2投影データに含まれる被検体サイノグラムのうち、復元の対象とするエネルギー帯の被検体サイノグラムを再構成して復元画像を生成する。
以上のステップS21〜S25の動作によって、画像処理部34bによる画像処理が実行される。
以上のように、生成部342aによって、正則化項および第1重みwjが加味された式(20)を用いて第2投影データが算出されることにより、第2投影データの復元スペクトルは、検出器13により検出された検出スペクトルと比較して、入射スペクトル(被検体スペクトル)においてエネルギー方向に急峻に変化している部分についても近似した形状に近づけることができる。
<変形例>
本変形例の画像処理部34cについて、第2の実施形態の画像処理部34bと相違する点を中心に説明する。本変形例の画像処理部34cは、算出部が、再構成部によって算出された物質密度を使用して、第1重みwjを算出する動作を中心に説明する。
図14は、第2の実施形態の変形例の画像処理部のブロック構成の一例を示す図である。図14を参照しながら、本変形例の画像処理部34cのブロック構成および動作について説明する。
図14に示すように、画像処理部34cは、記憶部341cと、生成部342cと、再構成部343cと、算出部344cと、を備えている。
記憶部341cは、上述の式(6)に示す行列Hおよび重みλを検出器応答データとして記憶する機能部である。
生成部342cは、データ収集部16から、被検体40のサイノグラムである被検体サイノグラムを第1投影データとして受信し、記憶部341cから検出器応答データを読み出し、第1投影データおよび検出器応答データに基づいて、サイノグラムの形式である第2投影データを生成する機能部である。生成部342cは、後述する算出部344cにより算出された第1重みwj、および、記憶部341cから読み出した検出器応答データ(行列H、重みλ)から、上述の式(6)、(16)、および(18)〜(20)によって、第1投影データに含まれる歪みを補正した第2投影データを算出して生成する。
算出部344cは、生成部342cが第2投影データを算出するために、上述の式(16)の第1重みwjを算出する機能部である。算出部344cは、例えば、X線管11から出射されるX線のスペクトルである出射スペクトルの情報、X線管11のターゲットの部材、X線管11のターゲットの角度、X線管11が備えるフィルタ、撮影時のX線管11の管電圧、使用した造影剤、および、想定される被検体40の組成のうちの少なくともいずれかの情報に基づいて、第2投影データの復元スペクトルがエネルギー方向に急峻に変換すると想定される程度を求めて第1重みwjを算出する。さらに、算出部344cは、後述する再構成部343cにより算出された、被検体40に含まれる物質の物質密度に基づいて、第2投影データの復元スペクトルがエネルギー方向に急峻に変換すると想定される程度を求めて第1重みwjを算出する。
再構成部343cは、生成部342cにより生成された第2投影データに含まれる被検体サイノグラムのうち、復元の対象とするエネルギー帯の被検体サイノグラムを再構成して復元画像を生成する機能部である。再構成部343cによる再構成の方法は、第1の実施形態と同様である。
再構成部343cは、まず、被検体40に含まれると予想される物質を特定し、その物質の種類の数と同数以上の種類のエネルギー帯の被検体サイノグラムを第2投影データから抽出して、抽出した被検体サイノグラムを再構成して、画素値が線減弱係数である複数の復元画像を生成する。次に、再構成部343cは、第1の実施形態の変形例2と同様の方法で、生成した複数の復元画像の画素毎に、特定した物質密度を算出する。そして、再構成部343cは、算出した物質密度の情報を、算出部344cに送る。
算出部344cは、例えば、再構成部343cから受け取った物質密度により、被検体40にヨードが含まれていると判定した場合、33[keV]の周辺のエネルギーで第1重みwjが小さくなるように算出する。そして、生成部342cは、上述したように、算出部344cにより算出された第1重みwjを使用して第2投影データを生成し、再構成部343cは、生成部342cにより改めて生成された第2投影データに含まれる被検体サイノグラムのうち、復元の対象とするエネルギー帯の被検体サイノグラムを再構成して復元画像を生成する。
なお、算出部344cは、再構成部343cから受け取った物質密度の高低に応じて第1重みwjを算出するものとしてもよい。例えば、物質密度が高いほど、第2投影データの復元スペクトルのエネルギー方向の変化は大きくなるので、第1重みwjを小さくするものとしてもよい。
また、第1重みwjは、被検体サイノグラムの画素毎に個別に算出してもよい。例えば、算出部344cは、再構成部343cから物資密度を密度画像の形式で受け取るものとし、密度画像においてヨードを含む画素を特定し、その画素をX線の透過経路に含む第1投影データの画素についてのみ、ヨードのK吸収端のエネルギーの周辺で第1重みwjが小さくなるように算出するようにしてもよい。
以上のように、算出部344cは、X線管11から出射されるX線のスペクトルである出射スペクトルの情報等だけでなく、再構成部343cにより算出された物質密度に基づいて、第1重みwjを算出する。したがって、第2投影データの復元スペクトルは、検出器13により検出された検出スペクトルと比較して、入射スペクトル(被検体スペクトル)においてエネルギー方向に急峻に変化している部分について、さらに近似した形状に近づけることができる。
(第3の実施形態)
本実施形態の画像処理部について、第2の実施形態の画像処理部34bと相違する点を中心に説明する。本実施形態の画像処理部の生成部342dは、スペクトルのエネルギー方向の滑らかさの異なる2つのスペクトルについての投影データを生成し、それらを合成する。以下、この2つのスペクトルについての投影データの生成および合成の動作の詳細を説明する。
図15は、第3の実施形態の画像処理部の生成部のブロック構成の一例を示す図である。図15を参照しながら、本実施形態の画像処理部の生成部342dのブロック構成について説明する。なお、本実施形態の画像処理部のブロック構成は、第2の実施形態の画像処理部34bのブロック構成と同様である。
図15に示すように、生成部342dは、第1補助生成部3421と、第2補助生成部3422と、合成部3423と、を備えている。
第1補助生成部3421は、算出部344により算出された第1重みwj、および、記憶部341bから読み出した検出器応答データ(行列H、重みλ)から、式(6)、(16)、および(18)〜(20)によって、第1投影データに含まれる歪みを補正した第3投影データを算出して生成する機能部である。
第2補助生成部3422は、算出部344により算出された第1重みwj、および、記憶部341bから読み出した検出器応答データ(行列H、重みλ)から、式(6)、(16)、および(18)〜(20)によって、第1投影データに含まれる歪みを補正した第4投影データを算出して生成する機能部である。このとき、第2補助生成部3422は、第4投影データについての復元スペクトル(第4スペクトル)の方が、第3投影データについての復元スペクトル(第3スペクトル)よりも、エネルギー方向の滑らかさが高まるように生成する。すなわち、第2補助生成部3422は、式(16)で用いる第1重み(式(16)ではwjと表記)を、第1補助生成部3421が同じく式(16)で用いる第1重みよりも大きくして、第4投影データを生成する。一例としては、第1補助生成部3421で用いる第1重みに1以上の所定値を乗じた値を、第2補助生成部3422で用いる第1重みとする方法がある。別の一例としては、第2補助生成部3422で用いる第1重みが、第1補助生成部3421で用いる第1重みよりも大きくなるように、予め定めた固定値にするという方法もある。
合成部3423は、第1補助生成部3421により生成された第3投影データと、第2補助生成部3422により生成された第4投影データとを合成して第2投影データを生成する機能部である。合成部3423は、算出部344から受け取る第1重みw
jが小さいほど、第3投影データに大きな重みを与える。具体的には、合成部3423は、以下の式(22)を用いて、第2投影データを生成する。
式(22)におけるeは、エネルギー帯の値を示すインデックスを示す。xpは、合成部3423によって生成される第2投影データを示すデータであり、具体的には、第2投影データを構成するそれぞれの被検体サイノグラムの画素pの画素値(フォトン数)で構成されるベクトルである。xp,3は、第1補助生成部3421によって生成される第3投影データを示すデータであり、具体的には、第3投影データを構成するそれぞれの被検体サイノグラムの画素pの画素値(フォトン数)で構成されるベクトルである。xp,4は、第2補助生成部3422によって生成される第4投影データを示すデータであり、具体的には、第4投影データを構成するそれぞれの被検体サイノグラムの画素pの画素値(フォトン数)で構成されるベクトルである。gは、0〜1の範囲の値をとる重みであり、第1重みwjが小さいほど、小さい値に設定される。例えば、第1重みwjを所定の定数で除した値を、1以下の値に丸めた値を重みgとすればよい。
図16は、第3の実施形態の画像処理部の動作の一例を示すフローチャートである。図16を参照しながら、本実施形態の画像処理部の画像処理動作について説明する。
<ステップS31>
画像処理部の生成部342dは、データ収集部16から、被検体40のサイノグラムである被検体サイノグラムを第1投影データとして受信して取得する。そして、ステップS32へ移行する。
<ステップS32>
生成部342dは、画像処理部の記憶部341bから検出器応答データを読み出して取得する。そして、ステップS33へ移行する。
<ステップS33>
算出部344は、生成部342bが第2投影データを算出するために利用する上述の式(16)の第1重みwjを算出する。まず、算出部344は、例えば、X線管11から出射されるX線のスペクトルである出射スペクトルの情報、X線管11のターゲットの部材、X線管11のターゲットの角度、X線管11が備えるフィルタ、撮影時のX線管11の管電圧、使用した造影剤、および、想定される被検体40の組成のうちの少なくともいずれかの情報を、システム制御部36等から取得する。そして、算出部344は、取得した情報に基づいて、第2投影データの復元スペクトルがエネルギー方向に急峻に変換すると想定される程度を求めて第1重みwjを算出する。算出部344は、算出した第1重みwjを、生成部342dに送る。
<ステップS34>
第1補助生成部3421は、算出部344により算出された第1重みwj、および、記憶部341bから読み出した検出器応答データ(行列H、重みλ)から、式(6)、(16)、および(18)〜(20)によって、第1投影データに含まれる歪みを補正した第3投影データを算出して生成する。第1補助生成部3421は、生成した第3投影データを合成部3423に送る。そして、ステップS35へ移行する。
<ステップS35>
第2補助生成部3422は、算出部344により算出された第1重みwj、および、記憶部341bから読み出した検出器応答データ(行列H、重みλ)から、式(6)、(16)、および(18)〜(20)によって、第1投影データに含まれる歪みを補正した第4投影データを算出して生成する機能部である。このとき、第2補助生成部3422は、第4投影データについての復元スペクトルの方が、第3投影データについての復元スペクトルよりも、エネルギー方向の滑らかさが高まるように生成する。すなわち、第2補助生成部3422は、式(16)で用いる第1重み(式(16)ではwjと表記)を、第1補助生成部3421が同じく式(16)で用いる第1重みよりも大きくして、第4投影データを生成する。第2補助生成部3422は、生成した第4投影データを合成部3423に送る。そして、ステップS36へ移行する。
<ステップS36>
合成部3423は、第1補助生成部3421により生成された第3投影データと、第2補助生成部3422により生成された第4投影データとを合成して第2投影データを生成する機能部である。合成部3423は、上述の式(22)を用いて、算出部344から受け取った第1重みwjが小さいほど、第3投影データに大きな重みを与えるようにして、第2投影データを生成する。
以上のステップS31〜S36の動作によって、本実施形態の画像処理部による画像処理が実行される。
以上のように、生成部342dは、復元スペクトルのエネルギー方向の滑らかさを抑えた第3投影データと、復元スペクトルのエネルギー方向の滑らかさを高めた第4投影データとを合成して、第2投影データを生成するものとしている。これによって、第2投影データの復元スペクトルは、検出器13により検出された検出スペクトルと比較して、入射スペクトル(被検体スペクトル)においてエネルギー方向に急峻に変化している部分についても近似した形状に近づけることができる。
上述の実施形態およびその変形例に係る画像処理装置(コンソール装置30)は、CPU(Central Processing Unit)等のマイクロプロセッサと、ROM(Read Only Memory)と、RAM(Random Access Memory)と、HDD(Hard Disk Drive)等の外部記憶装置と、ディスプレイ等の表示装置と、キーボードまたはマウス等の入力装置と、を備えており、通常のコンピュータを利用したハードウェア構成となっている。
上述の実施形態およびその変形例に係る画像処理装置(コンソール装置30)で実行されるプログラムは、インストール可能な形式または実行可能な形式のファイルでCD−ROM(Compact Disk Read Only Memory)、フレキシブルディスク(FD)、CD−R(Compact Disk Recordable)、DVD(Digital Versatile Disk)等のコンピュータで読み取り可能な記録媒体に記録してコンピュータプログラムプロダクトとして提供されるように構成してもよい。
また、上述の実施形態に係るおよびその変形例に係る画像処理装置(コンソール装置30)で実行されるプログラムを、インターネット等のネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることにより提供するように構成してもよい。また、上述の実施形態およびその変形例に係る画像処理装置(コンソール装置30)で実行されるプログラムをインターネット等のネットワーク経由で提供または配布するように構成してもよい。
また、上述の実施形態およびその変形例に係る画像処理装置(コンソール装置30)で実行されるプログラムは、ROM等に予め組み込まれて提供されるものとしてもよい。
上述の実施形態およびその変形例に係る画像処理装置(コンソール装置30)で実行されるプログラムは、コンピュータを上述した画像処理装置の各部(生成部342、342a〜342d、再構成部343、343c、算出部344、344c、第1補助生成部3421、第2補助生成部3422および合成部3423)として機能させ得る。このコンピュータは、CPUがコンピュータ読取可能な記憶媒体からプログラムを主記憶装置上に読み出して実行することができる。なお、上述の画像処理装置の各部の一部または全部は、ソフトウェアであるプログラムではなく、ハードウェア回路によって実現されてもよい。
本発明のいくつかの実施形態およびその変形例を説明したが、これらの実施形態およびその変形例は、例として提示したものであり、発明の範囲を限定することは意図していない。これらの新規な実施形態およびその変形例は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、および変更を行うことができる。これらの実施形態およびその変形は、発明の範囲および要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。