以下に、図面を参照しながら、本実施の形態に係るX線CT(Computed Tomography)装置、情報処理装置、および情報処理方法を詳細に説明する。なお、以下の実施形態では、同一の参照符号を付した部分は同様の動作を行うものとし、重複する説明を適宜省略する。
図1は、本実施の形態に係るX線CT装置1の一例を示す図である。
X線CT装置1は、放射線の一例であるX線を被検体Pに照射する。X線CT装置1は、例えば、被検体Pの断面画像を得るスペクトラルCT装置や、フォトンカウンティングCT装置などである。
X線CT装置1は、架台装置10と、寝台装置20と、情報処理装置30と、を備える。架台装置10および寝台装置20は、情報処理装置30にデータや信号授受可能に接続されている。
架台装置10は、被検体Pに対してX線を照射し、被検体Pを透過したX線のスペクトル情報を収集する。架台装置10は、照射制御回路11と、X線発生器12と、回転フレーム15と、検出器13と、収集回路14と、を備える。
X線発生器12は、X線を発生し、被検体Pへ照射する。X線発生器12は、X線管12Aと、ウェッジ12Bと、コリメータ12Cと、を備える。
X線管12Aは、照射制御回路11から供給される電圧によって、X線を発生する真空管である。そして、X線管12Aは、発生したX線を被検体Pへ照射する。例えば、X線管12Aは、円錐状または角錐状の広がりを有するビーム状のX線を発生させる。
ウェッジ12Bは、X線管12Aから照射されたX線のX線量を調節するフィルターである。コリメータ12Cは、ウェッジ12BによってX線量の調節されたX線の照射範囲を絞り込むスリットである。
このため、X線発生器12から被検体Pへ照射されるX線のスペクトルは、X線管12Aに供給される高電圧、電流、線源に用いるターゲット(例えば、タングステン等)の種類、ターゲットアングル、フィルターの種類(たとえば、ベリリウム等)や厚み、および、ウェッジ12Bの種類や厚みなどによって定まる。
回転フレーム15は、リング状の支持部材である。回転フレーム15は、X線発生器12と検出器13とを、被検体Pを挟んで対向するように支持する。なお、回転フレーム15の円中心には、被検体Pが位置する。回転フレーム15は、被検体Pを回転中心として回転可能に設けられている。このため、X線発生器12および検出器13は、対向配置された状態を維持したまま、被検体Pを回転中心として、回転可能に配置されている。
検出器13は、X線管12Aから照射され被検体Pを透過したX線のスペクトルを検出する。言い換えると、検出器13は、被検体Pを透過したX線の、エネルギーごとの光子カウント数を示すスペクトルを検出する。
なお、本実施の形態では、検出器13で検出されたスペクトルを、検出スペクトルと称して説明する。検出スペクトルの詳細は後述する。
検出器13は、回転フレーム15に沿って周方向に回転しながら、X線発生器12の位置ごとに、検出スペクトルを検出する。X線発生器12の位置とは、被検体Pに対するX線発生器12(詳細にはX線管12A)の相対位置を示す。X線発生器12の位置は、「ビュー」と称される場合がある。例えば、X線発生器12の位置は、回転フレーム15の周方向の1周360°における、予め定めた基準位置を0°としたときの角度で表す。例えば、架台装置10が0.5°ごとに検出スペクトルを検出すると仮定する。この場合、X線発生器12の位置は、0.5°ごとの角度で表すことができる。
検出器13は、複数の検出素子16を備える。検出素子16は、入射したX線に応じた信号を出力する。検出素子16は、例えば、テルル化カドミウム(CdTe)系の半導体素子である。
本実施の形態では、検出器13は、X線発生器12に対向する対向面における、回転フレーム15の周方向(図1中、矢印Q方向参照)、および、該周方向に直交する直交方向(図1中、矢印Z方向参照)に沿って、二次元配列されている。
本実施の形態では、検出素子16は、X線の光子が1つ入射するごとに、入射した光子のエネルギーおよび光子数の計測が可能なパルス状の信号を、収集回路14へ出力する。
なお、検出器13は、直接変換型、および、間接変換型の何れであってもよい。直接変換型の検出器13とは、検出素子16に入射した光子を、直接、電気的な信号に変換する検出器13である。間接変換型の検出器13とは、検出素子16のX線入射側に、シンチレータを配置した構成である。
照射制御回路11は、X線発生器12および回転フレーム15の動作を制御する。照射制御回路11は、高電圧発生機能11Aと、調整機能11Bと、架台駆動機能11Cと、を備える。
高電圧発生機能11A、調整機能11B、架台駆動機能11Cによって行われる各処理機能は、コンピュータによって実行可能なプログラムの形態で記憶回路へ記憶されている。照射制御回路11は、記憶回路から該各プログラムを読出し、実行することで、該各プログラムに対応する機能を実現するプロセッサである。換言すると、該各プログラムを読み出した状態の照射制御回路11は、図1の照射制御回路11内に示された各機能を有することとなる。
なお、図1には、単一の照射制御回路11によって、高電圧発生機能11A、調整機能11B、および架台駆動機能11Cで行われる処理機能が実現されるものとして説明する。しかし、複数の独立したプロセッサを組み合わせて照射制御回路11を構成し、各プロセッサがプログラムを実行することにより機能を実現するものとしてもよい。
高電圧発生機能11Aは、高電圧を発生し、発生した高電圧をX線管12Aへ供給する。調整機能11Bは、コリメータ12Cの開口度および位置を調整する。コリメータ12Cの開口度および位置の調整によって、X線管12Aから被検体Pに照射されるX線の照射範囲が調整される。例えば、調整機能11Bは、コリメータ12Cの開口度を調整することによって、X線の照射範囲、すなわち、X線のファン角およびコーン角を調整する。
架台駆動機能11Cは、回転フレーム15を回転駆動するように制御する。架台駆動機能11Cによる制御によって、回転フレーム15は、被検体Pを中心とした円軌道に沿って回転する。このため、X線発生器12および検出器13は、回転フレーム15によって、被検体Pを挟んで対向する位置関係を維持された状態で、回転フレーム15の回転に伴って被検体Pの周囲を回転する。
なお、架台装置10は、X線発生器12および検出器13を回転させる構成に限定されない。
例えば、架台装置10は、X線発生器12のみを回転させる構成であってもよい。この場合、検出器13を、回転フレーム15の周方向の1周分に渡って検出素子16を配列した構成とすればよい。また、回転フレーム15は、X線発生器12を支持する構成とすればよい。
収集回路14は、検出器13に設けられた複数の検出素子16の各々から受信した信号を用いて、検出素子16の各々に入射したX線の光子の数をカウントする。また、収集回路14は、複数の検出素子16の各々から受信した信号の示すパルス状の信号波形に基づいた演算処理を行う。この演算処理により、収集回路14は、複数の検出素子16の各々について、信号の出力を引き起こした光子のエネルギーを計測する。
これらの処理により、収集回路14は、X線のエネルギーに対する光子カウント数を示す検出スペクトルを、検出素子16の各々ごとに収集し、情報処理装置30へ出力する。このため、検出器13で検出した検出スペクトルとは、詳細には、検出素子16から出力された信号を収集回路14で収集することによって得られるスペクトルである。
ここで、上述したように、X線発生器12は、被検体Pの周囲を回転し、被検体Pに対する異なる複数の位置からX線を照射する。このため、本実施の形態では、収集回路14は、X線発生器12の位置ごとに、検出器13に設けられた複数の検出素子16の各々から信号を受信する。
収集回路14は、X線発生器12の位置および検出素子16ごとに、X線のエネルギーに対する光子カウント数を収集する。そして、収集回路14は、X線発生器12の位置および検出素子16ごとに、X線のエネルギーに対する光子カウント数を示す検出スペクトルを、情報処理装置30へ出力する。このとき、収集回路14は、X線発生器12の位置と、検出素子16の識別情報と、対応する検出スペクトルと、を含む収集情報を、情報処理装置30へ出力する。
なお、本実施の形態では、検出素子16の識別情報として、検出器13における検出素子16の位置を用いる場合を一例として説明する。このため、以下では、検出素子16の識別情報を、検出素子16の位置、と称して説明する場合がある。
寝台装置20は、被検体Pを載置する。寝台装置20は、寝台駆動装置21と、天板22とを、備える。
天板22は、被検体Pが載置されるベッド等の寝台である。天板22は、天板22に載置された被検体Pの体軸方向(図1中、矢印Z方向参照)に沿って移動可能に設けられている。寝台駆動装置21は、情報処理装置30の制御によって、天板22を被検体Pの体軸方向に沿って移動させる。なお、本実施の形態では、被検体Pの体軸方向と、回転フレーム15の回転軸方向と、が一致するものとして説明する。このため、天板22が被検体Pの体軸方向に沿って移動することで、天板22に載置された被検体Pは、回転フレーム15の内側に進入または該回転フレーム15の内側から外側へ搬送されることとなる。
情報処理装置30は、架台装置10および寝台装置20を制御する。
情報処理装置30は、インターフェース回路31と、ディスプレイ32と、入力回路33と、記憶回路34と、処理回路35と、を備える。
インターフェース回路31は、架台装置10および寝台装置20の各々と通信する。
入力回路33は、ユーザによる各種の操作指示を受付ける。入力回路33は、受付けた操作指示に応じた指示信号を、処理回路35へ出力する。入力回路33は、例えば、マウス、キーボード、ボタン、トラックボール、またはジョイスティックなどである。
ディスプレイ32は、各種画像を表示する。ディスプレイ32は、例えば、表示画像(詳細後述)や、CT画像(詳細後述)などを表示する。ディスプレイ32は、例えば、CRT(Cathode Ray Tube)ディスプレイ、LCD(Liquid Crystal Display:液晶ディスプレイ)、有機EL(Organic Electro−Luminescence)などである。
記憶回路34は、各種のデータを記憶する。記憶回路34は、例えば、HDD(Hard Disk Drive)、SSD(Solid State Drive)、RAM(Random Access Memory)、フラッシュメモリ等の半導体メモリ素子、光ディスクなどである。
処理回路35は、架台装置10および寝台装置20の制御や、架台装置10から取得した各種データを用いた情報処理を行う。
処理回路35は、スキャン制御機能35Aと、取得機能35Bと、演算機能35Cと、算出機能35Dと、更新機能35Eと、特定機能35Fと、生成機能35Gと、再構成機能35Hと、表示制御機能35Iと、受付機能35Jと、を備える。
スキャン制御機能35A、取得機能35B、演算機能35C、算出機能35D、更新機能35E、特定機能35F、生成機能35G、再構成機能35H、表示制御機能35I、および受付機能35Jの一部またはすべては、例えば、CPU(Central Processing Unit)などの電子回路にプログラムを実行させること、すなわち、ソフトウェアにより実現してもよい。
例えば、スキャン制御機能35A、取得機能35B、演算機能35C、算出機能35D、更新機能35E、特定機能35F、生成機能35G、再構成機能35H、表示制御機能35I、および受付機能35Jによって行われる各処理機能は、コンピュータによって実行可能なプログラムの形態で記憶回路34へ記憶されている。
情報処理装置30は、記憶回路34から該各プログラムを読出し、実行することで、該各プログラムに対応する機能を実現するプロセッサである。換言すると、該各プログラムを読み出した状態の処理回路35は、図1の処理回路35内に示された各機能を有することとなる。
なお、図1には、単一の処理回路35によって、スキャン制御機能35A、取得機能35B、演算機能35C、算出機能35D、更新機能35E、特定機能35F、生成機能35G、再構成機能35H、表示制御機能35I、および受付機能35Jで行われる処理機能が実現されるものとして説明する。しかし、複数の独立したプロセッサを組み合わせて処理回路35を構成し、各プロセッサがプログラムを実行することにより機能を実現するものとしてもよい。
また、スキャン制御機能35A、取得機能35B、演算機能35C、算出機能35D、更新機能35E、特定機能35F、生成機能35G、再構成機能35H、表示制御機能35I、および受付機能35Jの一部またはすべては、ハードウェアにより実現してもよいし、ソフトウェアおよびハードウェアを併用して実現してもよい。
該ハードウェアにより実現する場合、IC(Integrated Circuit)、ASIC(Application Specific Integrated Circuit)、FPGA(Field Programmable Gate Array)などの電子回路で実現すればよい。
スキャン制御機能35Aは、照射制御回路11、収集回路14、および、寝台駆動装置21を制御する。具体的には、スキャン制御機能35Aは、照射制御回路11を制御することによって、回転フレーム15を回転させ、X線管12AからX線を被検体Pへ照射させ、コリメータ12Cの開口度及び位置の調整を行う。この制御によって、架台装置10は、回転フレーム15を回転させながら、X線発生器12からX線を連続的または間欠的に被検体Pへ照射する。
例えば、スキャン制御機能35Aは、ヘリカルスキャン、または、ノンヘリカルスキャンを実行するように、架台装置10を制御する。
また、スキャン制御機能35Aは、収集回路14から収集情報を受付け、記憶回路34へ記憶する。上述したように、収集情報は、X線発生器12の位置と、検出素子16の位置と、対応する検出スペクトルと、を含む。
ここで、スペクトルについて説明する。図2(A)〜(D)は、スペクトルの一例の説明図である。
X線発生器12から照射されたX線は、被検体Pを透過して検出器13へ到り、検出器13で検出される。
以下では、X線管12Aから照射され、被検体Pへ到るまでのX線のスペクトルを、照射スペクトル50と称する。例えば、図2(A)に示すように、X線発生器12から検出素子16へX線a1が照射されたとする。この場合、X線a1における、領域XAのX線のスペクトルが、照射スペクトル50に相当する(図2(B)参照)。なお、領域XA上は空気領域であり、X線の減衰は無視できるものとする。
照射スペクトル50は、X線管12Aに供給される高電圧、電流、線源に用いるターゲット(例えば、タングステン等)の種類、および、ウェッジ12Bの種類や厚みなどによって定まる。このため、照射スペクトル50は、精密な実計測や、高精度なシミュレーションにより得られる。このため、本実施の形態では、記憶回路34が、予め照射スペクトル50を記憶する。
また、本実施の形態では、被検体Pを透過してから検出素子16へ到るまでのX線のスペクトルを、透過スペクトル52と称する。例えば、図2(A)に示すように、X線発生器12から検出素子16へ照射されたX線a1における、領域XBのX線のスペクトルが、透過スペクトル52に相当する(図2(C)参照)。
X線管12Aから被検体Pへ照射された照射スペクトル50は、被検体Pを透過することによって減衰する。この減衰は、被検体Pに含まれる物質に依存する。このため、被検体Pを透過したX線のスペクトルである透過スペクトル52は、被検体Pに含まれる物質に応じて減衰したスペクトル形状となる。すなわち、透過スペクトル52は、X線が被検体Pを透過した経路に存在する物質の組成と、被検体Pに照射された照射スペクトル50と、によって定まる。
また、本実施の形態では、検出器13で検出されたX線のスペクトルを、検出スペクトル54と称する。具体的には、検出器13で検出されたX線のスペクトルが、検出スペクトル54に相当する(図2(D)参照)。
ここで、理想的には、透過スペクトル52と、検出スペクトル54とは等しくなるはずであるが、実施に測定される検出スペクトル54は、現実の検出器13の検出器応答等を考慮すると、検出器13に入射される直前の透過スペクトル52とは、一般には異なったものになる。
なお、図2(B)〜図2(D)に示すスペクトル(照射スペクトル50、透過スペクトル52、検出スペクトル54)の形状は、一例であり、これらの形状に限定されない。
ここで、照射スペクトル50および検出スペクトル54の各々を、下記式(式(1)、式(2))の各々で定義する。
本実施の形態では、照射スペクトル50を、上記式(1)で定義する。式(1)中、cは、検出素子16の位置を示す。eは、X線のエネルギー(keV)を示す。
式(1)中、S0(c,e)は、照射スペクトル50を示す。具体的には、S0(c,e)は、位置cの検出素子16に照射されたX線の、エネルギーeの各々に対応する光子カウント数の群である。より具体的には、S0(c,e)は、被検体Pがなかったと仮定した場合に、位置cの検出素子16に、エネルギーeで照射されると考えられるX線の光子カウント数の推定値の群である。換言すると、照射スペクトル50は、例えば、最終的に位置cの検出素子16において検出される、X線のエネルギーe毎の光子カウント数の分布を示すX線スペクトルであって、X線発生器12と位置cの検出素子16とを結ぶ経路において、X線発生器12から照射されたX線が被検体Pに到達するまでのX線スペクトルである。なお、照射スペクトル50は、X線発生器12の位置に拘らず一定であることから、式(1)には、X線発生器12の位置を示すパラメータは含まれていない。但し、式(1)は、X線発生器12の位置を示すパラメータを更に含んだものであってもよい。なお、照射スペクトル50の例としては、上述の例に限られず、S0(e)のように、照射スペクトル50は、検出素子16の位置cを引数として持たなくてもよい。また、照射スペクトル50の例としては、S0のように、X線のエネルギーeについて積分された照射スペクトル50が用いられても良い。
また、本実施の形態では、検出スペクトル54を、上記式(2)で定義する。式(2)中、cは、検出素子16の位置を示す。eは、X線のエネルギー(keV)を示す。式(2)中、vは、X線発生器12の位置を示す。S1(c,v,e)は、検出スペクトル54を示す。詳細には、S1(c,v,e)は、X線発生器12が位置vにあるときに位置cの検出素子16で検出されたX線の、エネルギーeの各々に対応する光子カウント数の群である。
式(1)および式(2)中の、eについて具体的に説明する。本実施の形態では、照射スペクトル50および検出スペクトル54は、各々、予め定めたエネルギー幅のエネルギー帯(以下、単に、エネルギーと称する場合もある)ごとの光子カウント数を示す。
検出スペクトル54によって示される、光子カウント数の各々に対応するエネルギー帯のエネルギー幅を、第1のエネルギー幅E1と称する(図2(D)参照)。
また、照射スペクトル50によって示される、光子カウント数の各々に対応するエネルギー帯のエネルギー幅を、第2のエネルギー幅E2と称する(図2(B)参照)。
本実施の形態では、第2のエネルギー幅E2は、第1のエネルギー幅E1以下である。
すなわち、検出スペクトル54は、0keVから最大エネルギーに向かって、第1のエネルギー幅E1ごとのエネルギー帯の各々の光子カウント数を表したものである。なお、該最大エネルギーには、例えば、検出器13で検出可能な光子のエネルギーの最大値などを定めればよい。
また、照射スペクトル50は、0keVから上記最大エネルギーに向かって、第2のエネルギー幅E2ごとのエネルギー帯の各々の光子カウント数を表したものである。
なお、各エネルギー帯は、少なくとも一部が互いに重複していてもよいし、重複していなくてもよい。
第1のエネルギー幅E1および第2のエネルギー幅E2の値は限定されないが、20keV未満であることが好ましく、1keV以下であることが更に好ましい。
本実施の形態では、一例として、第1のエネルギー幅E1および第2のエネルギー幅E2が、1keVであるものとして説明する。
すなわち、本実施の形態では、照射スペクトル50および照射スペクトル50の各々は、1keVのエネルギー幅ごとに、各エネルギー幅のエネルギー帯に対応する光子カウント数を示す。
また、本実施の形態では、式(1)および式(2)中、eは、対応するエネルギー幅のエネルギー帯によって示されるエネルギーの上限値を示すものとして説明する。このため、本実施の形態では、e=1は、0keV以上1keV未満のエネルギー帯を示す。また、e=2は、1keV以上2keV以下のエネルギー帯を示す。同様に、e=100は、99keV以上100keV以下のエネルギー帯を示すものとして説明する。
次に、被検体Pを透過したX線の経路において生じるスペクトルの歪みや減衰について説明する。ここで、スペクトルの歪みとは、現実に検出される検出スペクトル54の、系が理想的である場合に期待されるスペクトルからのずれを意味する。ここで、系が理想的である場合に期待されるスペクトルとは、例えば照射スペクトル50と、被検体Pに含まれる物質の情報から、ランベルト・ベールの法則など単純な計算を用い、かつ検出器13の応答が理想的であると仮定した場合に期待されるスペクトルである。
スペクトルの歪みの第1の原因として、X線発生器12の特性に起因する照射スペクトル50の歪みが挙げられる。例えば、X線発生器12のX線管12aから照射されるX線の照射スペクトル50が、X線照射の焦点移動や、X線管12a側に設けられたコリメータの特性に起因して、歪む。
また、スペクトルの歪みの第2の原因として、X線発生器12から被検体Pまでの経路における照射スペクトル50の歪みが挙げられる。例えば、X線発生器12のX線管12aから照射されるX線の照射スペクトル50が、その照射経路において、例えばウェッジの、均一な照射スペクトル50を生成するためのCuフィルタのメカ精度のばらつきなどに起因して、歪む。
また、スペクトルの歪みの第3の原因として、X線が被検体Pを透過する際におけるビームハードニングに起因する透過スペクトル52の歪みが挙げられる。ここで、ビームハードニングとは、X線が物質を透過する際、物質のX線吸収率はエネルギー依存性があることから、波長の長くエネルギーの低いX線が、エネルギーの高いX線と比較して多く吸収されることから、X線の線質が変化し、その結果、物質質量減弱係数に線形性が失われることを言う。ビームハードニングにより、単純な方法を用いて算出されれたスペクトルからのずれが、透過スペクトル52に表れる。このようにして、被検体透過によるX線の減衰により、被検体Pを透過したX線の経路において、スペクトルの歪が生じる。かかる場合、被検体透過によるX線の減衰により、スペクトルの歪みが生じる。なお、ビームハードニングがエネルギースペクトル歪として精密に測定あるいは推定できれば、そのビームハードニングの程度(あるいはビームハードニングのみによるスペクトル歪)は、被検体のその透過経路の物質の組成を表すと言える。
また、スペクトルの歪みの第4の原因として、X線が被検体Pを透過する際における散乱に起因する透過スペクトル52の歪みが挙げられる。通常、X線管12Aから被検体Pへ照射されたX線は、照射された進行方向の検出器13において検出される。しかしながら、散乱が生じる場合、X線管12Aから被検体Pへ照射されたX線は、被検体P上で進行方向が変化し、当初照射された進行方向とは異なる角度に位置する検出器13において検出される。この結果、透過スペクトルに歪みが生じる。
また、スペクトルの歪みの第5の原因として、検出器13側のコリメータ(e.g. post-patient filter)等による透過スペクトル52の歪みが挙げられる。
また、スペクトルの歪みの第6の原因として、検出器13の応答特性が挙げられる。すなわち、検出器13で実際に検出される検出スペクトル54は、透過スペクトル52を更に歪ませた形となる。ここで、検出器13の応答特性とは、詳細には、検出器13に設けられた検出素子16の各々の応答特性である。検出素子16の応答特性とは、検出素子16に入射したX線の透過スペクトル52(あるいは検出スペクトル54)に歪みを生じさせる要因である。検出素子16の応答特性は、例えば、検出素子16に入射したX線のエネルギーごとの、エスケープ、蛍光、クロストーク、および散乱の発生確率や検出されるエネルギーのばらつき、などである。この場合、X線検出器13の応答特性により、被検体Pを透過したX線の経路においてスペクトルの歪みが生じる。
図1に戻り、そこで、本実施の形態のX線CT装置1では、情報処理装置30の処理回路35が、特有の機能を備える。以下、各機能について詳細を説明する。
以下、全体の処理の流れの概略を説明したのち、各処理の詳細について説明を行う。ここで、処理回路35が行う一連の処理の目的は、X線が被検体Pに含まれる弁別対象の物質を透過した距離(X線透過距離)を正確に算出することであるが、かかるX線透過距離の算出を検出スペクトル54から直接解法により求めるのは、前述のスペクトルの歪みが問題となり難しい場合もある。
そこで、実施形態においては、処理回路35は、演算機能35Cにより、X線が被検体Pに含まれる弁別対象の物質のX線透過距離の暫定的な推定値を示す推定距離と、被検体Pを透過したX線の経路において生じるスペクトルの歪みを示す情報とに基づいて所定のシミュレーションを行って、X線が被検体Pを透過した後のスペクトルとして推定される推定スペクトルを計算する。続いて、処理回路35は、算出機能35Dにより、検出スペクトル54の理論値である推定スペクトルを、検出器13で検出された、被検体を透過した後のスペクトルである検出スペクトル54の実測値と比較する。具体的には、処理回路35は、算出機能35Dにより、推定スペクトルと、検出スペクトル54との誤差を算出する。
ここで、処理回路35は、典型的には、例えば、2以上の推定距離の候補に対して、それぞれ推定スペクトルを算出し、算出された推定スペクトルを基に、検出スペクトル54との誤差を算出する。続いて、処理回路35は、特定機能35Fにより、算出された誤差に基づいて、最も真の値に近いと考えられるX線透過距離を特定する。例えば、処理回路35は、特定機能35Fにより、2以上の推定距離の候補に対して算出された誤差の中から、誤差が最小となるような推定距離を、物質のX線透過距離として特定する。
続いて、処理回路35は、更新機能35Eにより、特定したX線透過距離を基に、推定距離の値を更新し、次のイテレーションにおける引数とする。
このように、処理回路35は、特定機能35により、推定距離の更新、更新後の推定距離に基づいて推定スペクトルの演算を演算機能35Cに行わせる処理、誤差の算出を算出機能35Dにより行わせる処理、及びX線透過距離の特定の一連の処理を、所定の条件が満たされるまで繰り返す。かかるのち、処理回路35は、特定機能35Fにより、所定の条件が満たされた時の推定距離を、物質のX線透過距離として特定する。このようにして得られた物質のX線透過距離が、最終的に求めたい物質のX線透過距離となる。
なお、処理回路35は、例えば軟組織と造影剤等、2種類以上の物質に対して、上述の処理を行っても良い。かかる場合、処理回路35は、当該2種類以上の物質それぞれに対して、複数の推定距離の候補に対して推定スペクトルを算出し、算出した推定スペクトルを基に、誤差が最小となるような推定距離を、物質のX線透過距離として特定する。
次に、処理回路35の有する各機能及び、処理回路35が行う処理の手順について、より詳細に説明する。
受付機能35Jは、ユーザによる各種操作指示を入力回路33から受け付ける。
取得機能35Bは、検出スペクトル54と、推定距離の初期値と、を取得する。
検出スペクトル54は、上述したように、検出器13で検出されたスペクトルである。取得機能35Bは、スキャン制御機能35Aによって記憶回路34に記憶された収集情報を読取る。これによって、取得機能35Bは、X線発生器12の位置の各々において、複数の検出素子16の各々で検出された検出スペクトル54を取得する。
推定距離は、X線が被検体Pに含まれる弁別対象の物質を透過した距離(有効距離)の推定値を示す。具体的には、推定距離とは、X線発生器12から照射され被検体Pを透過して検出器13に入射したX線が、該X線の光路において、弁別対象の物質を透過した距離の合計値を推定した値である。
ここで、推定距離は、後述する更新機能35Eによって更新される。取得機能35Bが取得する推定距離の初期値は、更新機能35Eによって更新された推定距離ではなく、別途取得した推定距離である。
なお、弁別対象の物質が複数種類である場合、取得機能35Bは、複数種類の物質の各々の、推定距離の初期値を取得する。
取得機能35Bは、例えば、入力回路33から、弁別対象の物質の推定距離の初期値を取得する。なお、本実施の形態では、弁別対象の物質が複数種類である場合を説明する。このため、取得機能35Bは、入力回路33から、複数の物質の各々の推定距離の初期値を、入力回路33から取得する。
記憶回路34は、複数の物質の各々の種類を示す種類情報と、推定距離の初期値と、を対応づけて予め記憶してもよい。この場合、例えば、ユーザは、入力回路33を操作することで、弁別対象の物質の種類を入力する。受付機能35Jは、入力回路33から、弁別対象の物質の種類を示す種類情報の入力を受け付ける。そして、取得機能35Bは、受付けた種類情報に対応する推定距離の初期値を記憶回路34から読取る。これにより、取得機能35Bは、弁別対象の物質の各々の、推定距離の初期値を取得してもよい。
弁別対象の物質は、該物質を透過することでスペクトルの減衰が生じる物質であればよい。弁別対象の物質は、原子や分子であってもよいし、人体の特定の部位や組成(骨、筋、など)であってもよい。弁別対象の物質は、例えば、水、ヨード、カルシウム、ガドリニウム、筋肉、脂肪などである。
弁別対象の物質は、1種類であってもよいし、複数種類であってもよい。本実施の形態では、弁別対象の物質は、複数種類である場合を説明する。また、本実施の形態では、弁別対象の物質が、水、ヨード、カルシウム、の3種類である場合を一例として説明する。
弁別対象の物質の種類は、例えば、ユーザによる入力回路33の操作によって入力される。なお、記憶回路34に複数種類の物質の各々を示す種類情報を予め記憶してもよい。そして、ユーザによる入力回路33の操作指示により、該種類情報の中からユーザ所望の種類情報の物質を、弁別対象として選択させてもよい。また、取得機能35Bが、カルテや過去の検査情報、検査依頼情報、検査部位などの事前情報を記憶回路34などから取得し、該事前情報から弁別対象の物質を自動設定してもよい。
なお、詳細は後述するが、推定距離は、更新機能35Eによって、より最適な値となるように更新されていく。このため、推定距離の初期値は、被検体Pに含まれる物質の実際のX線透過距離により近い値であることが好ましい。このため、取得機能35Bは、弁別対象の物質の推定距離の初期値を、算出することによって取得することが好ましい。
取得機能35Bが初期値を算出することによって、後述する一連の処理の繰り返し数の削減や、後述する特定機能35Fが特定する物質のX線透過距離の精度向上や安定性の向上を図ることができる。
例えば、取得機能35Bは、公知の方法を用いて、弁別対象の物質の推定距離の初期値を算出することで、該初期値を取得してもよい。例えば、取得機能35Bは、検出スペクトル54と、弁別対象の物質の各々の平均吸収係数と、を用いて、特開平5−161633号公報に開示された方法を用いて、物質の厚さを算出してもよい。そして、取得機能35Bは、算出した物質の厚さを、弁別対象の物質の推定距離の初期値として用いてもよい。
なお、取得機能35Bは、例えば、照射スペクトル50および検出スペクトル54からX線CT画像を生成し、生成したX線CT画像から簡易的に初期値を算出してもよい。例えば、取得機能35Bは、照射スペクトル50および検出スペクトル54のそれぞれの、合計光子カウント数または実効エネルギーから、X線CT画像を生成する。このとき、照射スペクトル50および検出スペクトル54の各々の全エネルギー帯を利用して、X線CT画像を生成することが好ましい。これにより、検出器13の応答特性の影響がキャンセルされたり、高いSN比の情報を利用することができる。そして、取得機能35Bは、該X線の光路におけるCT値(線減弱係数)の総和を、水の平均線減弱係数で除した値を、弁別対象物質の一つである水の初期値とする。そして、取得機能35Bは、残りの弁別対象物質の初期値を、ゼロとする。このことにより、被検体が全て水だと仮定した場合に、簡易的に推定した初期値を得ることができる。
また、取得機能35Bは、X線CT画像におけるCT値を閾値処理することで、骨や造影剤などの高輝度領域のみを抽出してもよい。そして、取得機能35Bは、抽出した領域の該X線の光路におけるCT値の総和を、ヨードやカルシウムに割り当てることで簡易的に推定した、初期値を用いてもよい。なお、X線CT画像では、固定の閾値を用いた場合であっても、比較的高精度に空気領域、脂肪領域、筋肉・血液領域、骨・造影剤領域の4領域に分類できる。また、CT値ヒストグラム情報を利用して、動的に閾値を決定する方法なども多数提案されている。
なお、取得機能35Bが、特開平5−161633号公報に開示された方法を用いて、推定距離の初期値を算出する場合には、X線発生器12の複数の位置と、複数の検出素子16の位置と、を複数のグループに分類し、グループごとに、算出を行うことが好ましい。グループごとに算出を行うことにより、SN比の向上による安定性の向上や、算出数の削減による算出時間の短縮を図ることができる。
演算機能35Cは、推定スペクトルを演算する。推定スペクトルは、被検体Pを透過したX線の、検出器13により検出されるスペクトルを推定したものである。
演算機能35Cは、照射スペクトル50と、推定距離と、応答情報と、に基づいて、推定スペクトルを演算する。
応答情報は、検出器13の応答特性等によって生じるスペクトルの歪みを示す情報である。応答特性については、上述したため、ここでは説明を省略する。応答情報は、予め記憶回路34に記憶すればよい。例えば、記憶回路34は、検出素子16の識別情報と、応答情報と、を対応付けて予め記憶する。
照射スペクトル50は、予め記憶回路34に記憶すればよい。例えば、記憶回路34は、X線発生器12の位置および検出素子16の位置、に対応する照射スペクトル50を予め記憶する。なお、照射スペクトル50は、X線発生器12の位置に拘らず、一定である。このため、記憶回路34は、検出素子16の位置に対応する種類の照射スペクトル50のみを予め記憶してもよい。
演算機能35Cは、X線発生器12の位置の各々について、複数の検出素子16の各々ごとに、推定距離と、応答情報と、照射スペクトル50と、を用いて、推定スペクトルを演算する。
図3は、処理回路35が行う処理の概要の説明図である。
例えば、演算機能35Cが、図3(A)に示す照射スペクトル50を演算に用いたと仮定する。また、照射スペクトル50によって示されるX線が被検体Pを透過した透過スペクトル52が、図3(B)に示すスペクトルであったと仮定する。そして、検出素子16で実際に検出された検出スペクトル54が、図3(C)に示すスペクトルであったと仮定する。
この場合、演算機能35Cは、照射スペクトル50と、推定距離と、応答情報と、に基づいて、推定スペクトル56を演算する(図3(D)参照)。
詳細には、演算機能35Cは、照射スペクトル50を、弁別対象の物質の線減弱係数および推定距離に応じて減衰させ、且つ、応答情報によって示される歪みに応じて歪ませることによって、推定スペクトル56を演算する。
物質の線減弱係数とは、弁別対象の物質の、X線のエネルギーに対する線減弱係数を示す。X線エネルギーごとの物質の線減弱係数は、予め、記憶回路34に記憶すればよい。例えば、記憶回路34は、物質の種類を示す種類情報と、X線エネルギーと、線減弱係数と、を対応づけて予め記憶する。演算機能35Cは、推定スペクトル56の演算時に、記憶回路34から、弁別対象の物質の種類情報に対応する線減弱係数を読取ることで、推定スペクトル56の演算を行えばよい。
まず、スペクトルの歪みが生じる原因が、前述の第6の原因である場合、すなわち、検出器13の応答特性に基づくスペクトルの歪みを考慮した推定スペクトル56を演算する場合について簡単に説明する。かかる場合、処理回路35は、演算機能35Cにより、スペクトルの歪みを表す情報である検出器13の応答特性と、推定距離に基づいて、推定スペクトル56を演算する。
具体的には、演算機能35Cは、下記式(3)により、推定スペクトル56を演算する。
式(3)中、cは、式(1)および式(2)と同様である。また、式(3)中、eは、式(1)と同様である。また、式(3)中、vは、式(2)と同様である。
また、式(3)中、mは、弁別対象の物質に種類ごとに通し番号を付与したときの該通し番号を示す。NMは、弁別対象の物質の数(種類数)を示す。上述したように、本実施の形態では、弁別対象の物質が、水、ヨード、カルシウム、の3種類である場合を一例として説明する。このため、m=1が水、m=2がヨード、m=3がカルシウムをそれぞれ表し、NM=3であるものとして説明する。
式(3)中、μmは、弁別対象の物質mの線減弱係数である。lmは、弁別対象の物質mの推定距離を示す。Fc,eは、検出器13の応答情報である。詳細には、Fc,eは、位置cの検出素子16の、X線のエネルギーeに対応する応答情報を示す。
S2(c,v,e)は、推定スペクトル56を示す。詳細には、S2(c,v,e)は、X線発生器12が位置vにあるときに位置cの検出素子16で検出されたX線の、エネルギーeに対応する光子カウント数の推定値の群である。
すなわち、式(3)の右辺における、Fc,e以外の部分が、So(c,e)によって表される照射スペクトル50を、弁別対象の物質の各々の推定距離lmおよび線減弱係数により減衰させたスペクトルを示す。そして、このスペクトルを応答情報Fc,eによって示される歪みに応じて歪ませることによって、S2(c,v,e)によって示される推定スペクトル56が算出される。
ここで、上述したように、式(3)中、eは、式(1)と同様である。すなわち、推定スペクトル56の演算時に用いる、エネルギーeのエネルギー幅は、照射スペクトル50の上記第2のエネルギー幅E2(図2(B)参照)と一致する。
このため、上述したように、本実施の形態では、式(3)におけるe=1は、0keV以上1keV未満のエネルギー帯を示す。また、e=2は、1keV以上2keV以下のエネルギー帯を示す。同様に、e=100は、99keV以上100keV以下のエネルギー帯を示すものとして説明する。
なお、演算機能35Cは、線減弱係数を記憶回路34から読取り、上記演算に用いればよい。このため、記憶回路34は、推定スペクトル56における各光子カウント数に対応するエネルギー帯の各々の、線減弱係数を予め記憶すればよい。例えば、記憶回路34は、各エネルギー帯の上限値と下限値との中央値の各々に対応する線減弱係数を、弁別対象の物質ごとに予め記憶すればよい。エネルギー帯の上限値と下限値との中央値は、例えば、エネルギー帯が99keV以上100keV以下の範囲である場合、99.5keVである。
なお、上述したように、本実施の形態では、第1のエネルギー幅E1および第2のエネルギー幅E2が同じエネルギー幅であるものとして説明している。このため、本実施の形態では、演算機能35Cは、上記式(3)により、推定スペクトル56を演算する場合を説明する。しかし、第1のエネルギー幅E1が第2のエネルギー幅E2より大きい場合、上記式(3)の右辺について、第1のエネルギー幅E1に対応するように、エネルギー帯の正規化(平均操作)を行えばよい。なお、この正規化には、公知の方法を用いればよい。
なお、応答情報は、単位時間あたりに検出器13に入射する光子の数によって変化する場合がある。例えば、光子が連続して検出器13に入射した場合、検出器13の反応が不十分になり、見かけ上、実際より低エネルギー(または高エネルギー)の光子の数としてカウントされる場合がある。このため、演算機能35Cは、検出スペクトル54によって示される光子カウント数に応じた応答情報を用いて、推定スペクトル56を演算してもよい。
この場合、記憶回路34は、検出素子16の位置と、光子カウント数と、の各々に対応する応答情報を予め記憶すればよい。
次に、スペクトルの歪みが生じる原因が、前述の第1〜第5の原因である場合に、スペクトルの歪みを考慮した推定スペクトル56を演算する場合について簡単に説明する。
スペクトルの歪みが生じる原因が、前述の第1の原因(X線照射の焦点移動、管球側のコリメータの特性等)または第2の原因(ウェッジ等)など、管球側の原因である場合には、処理回路35は、式(3)において、s0(c,e)を、F1(sO(c,e))に置き換えた式を用いて、式(3)の右辺を評価することにより、演算機能35Cにより、推定スペクトル56を演算する。ここで、F1は、前述の第1の原因または第2の原因により生じる照射スペクトル50の歪みを表す関数である。かかる関数の具体的な表式は、例えば所定のテーブルの形で、記憶回路34に保持される。なお、前述の第6の原因による補正を行わない場合は、Fc,e=1である。
また、スペクトルの歪みが生じる原因が、前述の第5の原因(検出器側のコリメータ等)、すなわち検出器側の原因である場合には、処理回路35は、式(3)におけるFc,eの表式の中に、これらの効果をとりこみ、式(3)の右辺を評価することにより、演算機能35Cにより、推定スペクトル56を演算する。かかる関数の具体的な表式は、例えば所定のテーブルの形で、記憶回路34に保持される。
また、スペクトルの歪みが生じる原因が、前述の第3の原因(ビームハードニング)である場合には、処理回路35は、例えばモンテカルロシミレーションを用いて、ビームハードニングによるスペクトル変化を示すテーブルを生成することができる。かかるテーブルにより、被検体の物質組成に起因する、X線減弱のエネルギー依存を知ることができ、これを基に、物質組成の推定が可能になる。処理回路35は、生成したテーブルを記憶回路34に保存する。処理回路35は、演算機能35Cにより、被検体Pの撮影に係る処理の際、記憶回路34に保存されたかかるテーブルと、推定距離に基づいて、推定スペクトル56を演算する。かかる場合、前述のテーブル(被検体透過によるX線の減衰を表す情報)が、スペクトルの歪みを表す情報の例となる。なお、かかる場合、実施形態は、モンテカルロシミュレーションをすることを要せず、より簡便な方法でかかるテーブルを作成してもよい。
また、スペクトルの歪みが生じる原因が、前述の第4の原因(散乱)である場合には、X線管12Aは、例えばファントムなどの既知の物質に対して、既知の照射スペクトル50のX線を照射する。検出器13は、検出スペクトル54を検出する。かかる検出スペクトル54には、散乱の効果が含まれているので、処理回路35は、ファントムを構成する物質の情報と、照射スペクトル50と、検出スペクトル54とに基づいて、散乱の効果を表す散乱情報テーブルを生成することができる。処理回路35は、生成した散乱情報テーブルを記憶回路34に保存する。処理回路35は、演算機能35Cにより、被検体Pの撮影に係る処理の際、記憶回路34に保存されたかかる散乱情報テーブルと、推定距離に基づいて、推定スペクトル56を演算する。かかる場合、前述の散乱情報テーブル(X線の散乱特性)が、スペクトルの歪みを表す情報の例となる。
次に、算出機能35Dについて説明する。算出機能35Dは、推定スペクトル56と検出スペクトル54との誤差を算出する。
具体的には、算出機能35Dは、X線管12Aの各位置の各々について、推定スペクトル56と検出スペクトル54との誤差を、検出素子16の各々ごとに算出する。算出機能35Dは、X線発生器12の位置および検出素子16の位置が同じ推定スペクトル56と検出スペクトル54との組合せごとに、誤差を算出する。
なお、後述する更新機能35Eによる推定距離の更新と、演算機能35Cによる推定スペクトル56の演算と、算出機能35Dによる誤差の算出と、の一連の処理は、後述する特定機能35Fの制御によって繰り返し実行される(詳細後述)。
このため、算出機能35Dは、演算機能35Cで直前に演算された推定スペクトル56を、誤差の算出対象として用いる。すなわち、算出機能35Dは、演算機能35Cで新たに推定スペクトル56が演算される度に、新たに演算された推定スペクトル56を用いて、検出スペクトル54との誤差を算出する。
例えば、算出機能35Dは、検出スペクトル54と推定スペクトル56との、エネルギーの各々に対する光子カウント数の差に基づいて、誤差を算出する。より具体的には、算出機能35Dは、検出スペクトル54と推定スペクトル56との、エネルギーの各々に対する光子カウント数の差の積算値(例えば二乗平均)を、誤差として算出する。
詳細には、算出機能35Dは、検出スペクトル54と推定スペクトル56との、第2のエネルギー幅E2のエネルギー帯の各々に対応する、光子カウント数の差の積算値を、誤差として算出する。
具体的には、例えば、算出機能35Dは、下記式(4)を用いて、誤差を算出する。
式(4)中、dは、誤差を示す。また、式(4)中、S1(c,v,e)は、上記式(2)と同様である。式(4)中、S2(c,v,e)は、上記式(3)と同様である。また、式(4)中、c、e、vの各々は、式(3)と同様である。
また、式(4)中、NEは、エネルギーeの上限値である。本実施の形態では、一例として、NE=120であるものとして説明する。
なお、算出機能35Dは、下記式(5)を用いて、誤差を算出してもよい。
式(5)中、d、S1(c,v,e)、S2(c,v,e)、c、e、および、vの各々は、式(4)と同様である。式(5)中、Eminは、検出器13で検出可能な光子のエネルギーの最小値より大きい値を示す。Emaxは、検出器13で検出可能な光子のエネルギーの最大値より小さい値を示す。
すなわち、式(5)に示すように、算出機能35Dは、誤差の算出時に用いるエネルギーの範囲を、Emin以上Emax以下の任意のエネルギー範囲に限定してもよい。
また、算出機能35Dは、下記式(6)、式(7)、および式(8)の何れかを用いて、誤差を算出してもよい。
式(6)〜式(8)中、d、S1(c,v,e)、S2(c,v,e)、c、e、v、Emin、および、Emaxの各々は、上記式(5)と同様である。また、関数gはフィルタ関数として用いる所定の関数であり、演算子"*"は畳み込み演算を表す。
式(6)中、a(e)は、エネルギーeに対応する重み付け係数を示す。重み付け係数には、エネルギー毎に任意の値を予め定めればよい。
換言すると、算出機能35Dは、検出スペクトル54と推定スペクトル56との、エネルギーの各々に対する光子カウント数の差と、弁別対象の物質の特性に基づいたエネルギーに応じた重み付け係数とに基づいて、推定スペクトルと、検出スペクトルとの誤差を算出する。一例として、式(6)のように、算出機能35Dは、検出スペクトル54と推定スペクトル56との、エネルギーの各々に対する光子カウント数の差の二乗と、弁別対象の物質の特性に基づいたエネルギーに応じた重み係数との線形和に基づいて、推定スペクトルと、検出スペクトルとの誤差を算出する。
例えば、誤差の算出に重要と想定されるエネルギー帯を予め設定する。そして、算出機能35Dは、該エネルギー帯のエネルギーに、該エネルギー帯以外のエネルギーに比べて高い重み付け係数を予め設定すればよい。
一例として、算出機能35Dは、所定のエネルギーとのエネルギー差が第1の閾値以下のエネルギーにおける重み付け係数が、その他のエネルギーにおける重み付け係数より小さな重み付け係数となるような重み付け係数に基づいて、誤差を算出する。ここで、例えば所定のエネルギーとは、K吸収端に対応するエネルギーである。かかる重み付け係数を採用する理由は、例えば、以下の通りである。すなわち、K吸収端に近いエネルギー領域においては、検出スペクトル54に誤差が生じやすいので、これらのエネルギー領域の重み付け係数を小さくするのが望ましい。従って、算出機能35Dは、K吸収端付近のエネルギー領域の重み付け係数を小さくする。
また、別の例として、算出機能35Dは、所定のエネルギーとのエネルギー差が第1の閾値より大きく第2の閾値以下のエネルギーにおける重み付け係数が、その他のエネルギーにおける重み付け係数より大きな重み付け係数となるような重み付け係数に基づいて、誤差を算出する。ここで、例えば所定のエネルギーとは、K吸収端に対応するエネルギーである。かかる重み係数を採用する理由は、例えば、以下の通りである。すなわち、前述のK吸収端に近いエネルギー領域を除外するとすれば、K吸収端に近いエネルギーほど、X線透過距離の推定を精度よく行うことができる。従って、算出機能35Dは、K吸収端とのエネルギー差が、第1の閾値より大きく、かつ第1の閾値より大きい第2の閾値以下のエネルギーとなるようなエネルギー領域の重み付け係数を大きくする。換言すると、Kエッジ近傍のスペクトル情報は、物質弁別に有効なので大きな重み付け係数を採用することが望ましいが、Kエッジ最近傍(Kエッジとのエネルギー差が第1の閾値以下)においては、誤差が大きくなる可能性があるため、小さな重み付け係数を採用する。
また、第1の閾値は、「0」であってもよい。かかる場合、算出機能35Dは、所定のエネルギーとのエネルギー差が第2の閾値以下のエネルギーにおける重み係数が、その他のエネルギーにおける重み付け係数より大きな重み付け係数となるような重み付け係数に基づいて、誤差を算出する。
また、重み付け係数には、弁別対象の物質の特性に基づいたエネルギーに応じた値を定めてもよい。弁別対象の物質の特性は、例えば、物質のX線に対する減衰特性、および、K吸収端のエネルギー、の少なくとも1つを含む。
例えば、算出機能35Dは、弁別対象の複数種類の物質の内、少なくとも1種類の物質のX線に対する減衰特性によって示される、減衰率が閾値以上または閾値未満のエネルギー帯の領域を予め設定する。この閾値には、任意の値を設定すればよい。また、この閾値は、ユーザによる入力回路33の操作によって、適宜更新可能としてもよい。そして、算出機能35Dは、該エネルギー帯のエネルギーに、該エネルギー帯以外のエネルギーに比べて高い重み付け係数を予め設定すればよい。
また、例えば、弁別対象の物質のK吸収端(例えば、ヨードの場合、33keV付近)を含むエネルギー帯に、該エネルギー帯以外のエネルギーに比べて高い重み付け係数を設定してもよい。
すなわち、算出機能35Dは、検出スペクトル54と推定スペクトル56との、エネルギーの各々に対する光子カウント数の差の各々に、物質の特性に基づいたエネルギーに応じた重み付け係数を付与した後の積算値を、誤差として算出してもよい。
また、式(7)に示すように、算出機能35Dは、検出スペクトル54に含まれるノイズの影響などを考慮するため、関数gを更に用いた畳み込み演算を行うことで、誤差を算出してもよい。
なお、関数gは、例えば、平均値フィルターやガウシアンフィルターなどの線形フィルターや、εフィルターなどの非線形フィルターである。各フィルターのフィルタサイズ(範囲)やεなどのフィルタパラメータは、例えば、ヨードのk吸収端である33keV付近のエネルギーをぼかさないようにするなど、エネルギーによって切り替えるとより好適である。
すなわち、算出機能35Dは、式(7)における関数gを用いて、誤差算出対象の検出スペクトル54と推定スペクトル56における、33keV付近のスペクトルがぼけないように、フィルター範囲を調整しても良い。
式(8)中、bは、連続するエネルギーeを予め定めた数(2以上)ごとの複数の群に分類したときの、各群の通し番号である。各群の通し番号とは、例えば、各群をエネルギーeの小さい方から大きくなる方向へ向かって、昇順に数値を割当てたものを用いればよい。また、式(8)中、NBは、該群の数を示す。
また、式(8)中、E0およびE1は、各々、群bに属するエネルギーeの下限値および上限値の各々を示す。例えば、E0(1)=1、E1(1)=20、E0(2)=21、E1(2)=40、・・・であり、NB=6である。
すなわち、式(8)に示すように、算出機能35Dは、推定スペクトル56および検出スペクトル54の各々における、任意のエネルギー範囲(E0(b)以上E1(b)以下の範囲)の光子カウント数の総和の差を、誤差として算出してもよい。
更新機能35Eは、推定距離を更新する。詳細には、更新機能35Eは、算出機能35Dが前回の誤差算出時に用いた推定スペクトル56に対応する推定距離を、更新する。前回の誤差算出時に用いた推定スペクトル56に対応する推定距離とは、該推定スペクトル56の演算時に演算機能35Cで用いた推定距離である。なお、以下では、前回の誤差算出時に用いた推定スペクトル56に対応する推定距離を、前回用いた推定距離、と称して説明する場合がある。
更新機能35Eによる、前回用いた推定距離の更新方法については、後述する。
なお、更新機能35Eによる更新とは、記憶回路34に既に記憶されている推定距離を残したまま、該推定距離の値を変更した後の推定距離を更に記憶回路34へ記憶することを意味する。このため、更新機能35Eが推定距離を更新する度に、更新後の新たな推定距離が記憶回路34へ順次記憶されることとなる。
特定機能35Fは、算出機能35Dで算出された誤差に基づいて、物質のX線透過距離を特定する。
例えば、特定機能35Fは、算出機能35Dで算出された誤差が閾値以下となる推定距離を、物質のX線透過距離として特定する。
本実施の形態では、特定機能35Fは、推定距離の更新、更新後の推定距離に基づいた推定スペクトル56の演算、および誤差の算出、の一連の処理を繰り返すように、更新機能35E、演算機能35C、および算出機能35Dを制御する。
図3に示すように、特定機能35Fによる制御によって、上記一連の処理毎に、推定距離が更新される(図3(E)参照)。なお、図3(E)中、L1は、弁別対象の物質である水の推定距離の一例を示し、L2はヨードの推定距離の一例を示し、L3は、カルシウムの推定距離の一例を示す。
すなわち、特定機能35Fは、推定距離(図3(E)参照)の更新、照射スペクトル50(図3(A)参照)と更新後の推定距離とに基づいた推定スペクトル56(図3(D)参照)の演算、および、推定スペクトル56と検出スペクトル54(図3(C)参照)との誤差の算出、の一連の処理を、繰返し実行する。
図1に戻り、そして、特定機能35Fは、閾値以下の誤差の算出された一連の処理で用いた推定距離を、物質のX線透過距離として特定する。
すなわち、特定機能35Fは、該閾値以下の誤差が算出されるまで、推定距離の更新、更新後の推定距離に基づいた推定スペクトル56の演算、および誤差の算出、の一連の処理を繰り返すように、更新機能35E、演算機能35C、および算出機能35Dを制御する。
なお、特定機能35Fが用いる該閾値には、誤差や推定距離が最適化されたか否かの判別に用いることの可能な値を定めればよい。例えば、特定機能35Fが用いる該閾値には、一連の処理毎に算出された誤差の内の最小値を定めればよい。この場合、特定機能35Fは、一連の処理毎に算出された誤差の内の最小値を、閾値以下の誤差として特定する。そして、特定機能35Fは、該閾値以下の誤差の算出された一連の処理で用いた推定距離を、物質のX線透過距離として特定すればよい。すなわち、特定機能35Fは、誤差が最小となる推定距離を、物質のX線透過距離として特定してもよい。
なお、更新機能35Eは、上記一連の処理毎に、前回算出された誤差より小さい誤差が算出されるように、前回用いた推定距離を更新することが好ましい。更新機能35Eによる推定距離の更新方法には、例えば、貪欲法、最急降下法、および共役勾配法などの、公知の最適化方法を用いればよい。
なお、最適化を行うにあたって、典型的には、処理回路35は、算出機能35Dにより、2以上の推定距離の候補に対して誤差を算出する。
かかる場合、処理回路35は、特定機能35Fにより、例えば貪欲法を用いる場合、前述の2以上の推定距離の候補に対して計算された誤差のうち、誤差が最小値になるような推定距離の候補が、次のイタレーションステップにおける推定距離となる。
また、最急降下法など、イタレーションステップの更新に、微分係数の算出が必要となる場合、処理回路35は、特定機能35Fにより、前述の2以上の推定距離の候補に対して計算された誤差に対して例えば差分法を用いて微分係数を算出し、算出した微分係数を用いて、次のイタレーションステップにおける推定距離を算出する。
このように、特定機能35Fによる制御によって、上記一連の処理毎に、推定距離が更新される。このため、算出機能35Dは、誤差の算出時に用いる重み付け係数(上記式(6)および式(7)におけるa(e)参照)やエネルギー範囲であるEmin、Emax、E0(b)、E1(b)、フィルタ関数gを、更新された推定距離に応じて、上記一連の処理毎に変更してもよい。
例えば、弁別対象の物質の各々の更新後の推定距離の内、ヨードの推定距離が最も大きかったとする。この場合、検出スペクトル54と次回の一連の処理時に演算される推定スペクトル56における、ヨードのk吸収端である33keV付近の光子カウント数の差が、今回の一連の処理時より小さくなることが好ましい。
そこで、この場合、算出機能35Dは、上記式(6)〜式(7)の各々に示される重み付け係数a(e)について、33keV付近のエネルギー帯に対応する重み付け係数a(e)を、該エネルギー帯以外のエネルギーに比べて大きくすることが好ましい。
また、算出機能35Dは、上記式(5)でEminとEmaxによって示される任意のエネルギー範囲を、33keV付近に限定してもよい。
このように、算出機能35Dは、更新された推定距離に応じて重み付け係数を更新することで、例えば、33keV付近の光子カウント数の差が誤差に与える影響を大きくすることができる。
また、算出機能35Dは、式(8)における、E0およびE1によって表される、群bに属するエネルギーeの下限値および上限値によって規定されるエネルギー帯内に、33keVが含まれるように、設計することが好ましい。
一方、弁別対象の物質の各々の更新後の推定距離の内、ヨードの推定距離が最も小さかったと仮定する。この場合、算出機能35Dは、上記式(6)〜式(7)の各々に示される重み付け係数a(e)について、33keV付近のエネルギー帯に対応する重み付け係数a(e)を、該エネルギー帯以外のエネルギーに比べて小さくすることが好ましい。
この場合、算出機能35Dは、更新された推定距離に応じて重み付け係数を変更することで、例えば、33keV付近の光子カウント数の差が誤差に与える影響を小さくすることができる。よって、この場合、透過する光子の数が少なく、ノイズの影響を受けやすい33keV付近の光子カウント数の差が、次回算出される誤差に与える影響を小さくすることができる。すなわち、この場合、ノイズの影響による、誤差算出精度の低下を抑制することができる。
次に、再構成機能35Hについて説明する。再構成機能35Hは、取得機能35Bで取得した、X線発生器12の位置ごとに、検出素子16の各々で検出された検出スペクトル54を用いて、公知の再構成方法により、CT画像を再構成する。再構成方法は、例えば、逆投影処理である。また、逆投影処理は、例えば、FBP(Filtered Back Projection)法である。
生成機能35Gは、特定機能35Fで特定されたX線透過距離に応じた表示画像を生成する。
表示画像は、特定機能35Fで特定された、弁別対象の物質の各々のX線透過距離、該X線透過距離から算出した物質の密度、弁別対象の物質の種類、および、弁別対象の物質の原子番号の少なくとも1つを更に含む画像である。
例えば、生成機能35Gは、X線発生器12の位置の各々に対応する、検出器13に設けられた検出素子16ごとに特定された、弁別対象の物質(例えば、水、ヨード、カルシウム)の各々のX線透過距離を用いて、被検体Pに含まれる、弁別対象の物質の各々の、被検体Pにおける密度分布を算出する。密度分布の算出には、公知の方法を用いればよい。
そして、生成機能35Gは、例えば、水の密度分布をGチャネル、ヨードの密度分布をRチャネル、カルシウムの密度分布をBチャネルに割り当てたカラー画像を、表示画像として生成する。
また、生成機能35Gは、弁別対象の物質の各々の被検体Pにおける密度を、密度の高さに応じた色および濃度で示すカラー画像を、表示画像として生成してもよい。
また、生成機能35Gは、被検体Pにおける物質の密度を示す画像から、単色X線のCT画像や高エネルギー帯のみのX線CT画像を生成し、表示画像としてもよい。
また、生成機能35Gは、上記カラー画像と、再構成機能35Hで生成されたCT画像と、を合成した合成画像を、表示画像として生成してもよい。この場合、生成機能35Gは、CT画像と、カラー画像と、をαブレンディングすることによって、合成画像(表示画像)を生成すればよい。
表示制御機能35Iは、各種画像をディスプレイ32へ表示する制御を行う。本実施の形態では、表示制御機能35Iは、生成機能35Gで生成された表示画像や、再構成機能35Hで生成されたCT画像を、ディスプレイ32へ表示する制御を行う。
次に、処理回路35で実行する弁別処理の手順を説明する。図4は、処理回路35で実行する弁別処理の手順の一例を示すフローチャートである。
まず、取得機能35Bが、検出スペクトル54と、推定距離の初期値と、を取得する(ステップS100)。
次に、演算機能35Cが、照射スペクトル50と、ステップS100で取得した推定距離の初期値と、応答情報と、に基づいて、推定スペクトル56を演算する(ステップS102)。
次に、算出機能35Dが、ステップS102で算出された推定スペクトル56と、ステップS100で取得した検出スペクトル54と、の誤差を算出する(ステップS104)。そして、算出機能35Dは、ステップS104で算出した誤差を、ステップS100で取得した推定距離に対応づけて、記憶回路34へ記憶する(ステップS106)。
次に、更新機能35Eが、推定距離を更新する(ステップS108)。更新機能35Eは、記憶回路34に記憶されている最新の推定距離を更新する。そして、更新機能35Eは、更新後の推定距離を、記憶回路34へ新たに記憶する(ステップS110)。
次に、演算機能35Cが、照射スペクトル50と、ステップS108で更新された更新後の推定距離と、応答情報と、に基づいて、推定スペクトル56を演算する(ステップS112)。
次に、算出機能35Dが、ステップS112で算出された推定スペクトル56と、ステップS100で取得した検出スペクトル54と、の誤差を算出する(ステップS114)。そして、算出機能35Dは、ステップS114で算出した誤差を、ステップS108で更新した推定距離に対応づけて、記憶回路34へ記憶する(ステップS116)。
このため、記憶回路34には、ステップS108〜ステップS116の一連の処理ごとに、検出スペクトル54と推定スペクトル56との誤差と、弁別対象の物質の各々の推定距離と、が対応づけて順次記憶されていく。
次に、特定機能35Fが、処理を終了するか否かを判断する(ステップS118)。特定機能35Fは、例えば、ステップS108〜ステップS116の一連の処理を、予め定めた回数繰り返したか否かを判別することで、ステップS118の判断を行う。
なお、特定機能35Fは、処理終了を示す指示信号を入力回路33から受付機能35Jで受け付けたか否かを判別することで、ステップS118の判断を行ってもよい。また、特定機能35Fは、ステップS114で算出された誤差が、特定機能35Fが用いる上述した閾値以下の値であるか否かを判別することで、ステップS118の判断を行ってもよい。
また、特定機能35Fは、ステップS114で算出された誤差の、前回算出した誤差からの変動が予め定めた変動値以下となったか否かを判別することで、ステップS118の判断を行ってもよい。
ステップS118で否定判断すると(ステップS118:No)、ステップS108へ戻る。すなわち、特定機能35Fは、ステップS108〜ステップS118:Noの処理を繰り返すように制御する。すなわち、特定機能35Fは、推定距離の更新、更新後の推定距離に基づいた推定スペクトル56の演算、および、該推定スペクトル56を用いた誤差の算出、の一連の処理を繰り返すように、更新機能35E、演算機能35C、および算出機能35Dを制御する。
一方、ステップS118で肯定判断すると(ステップS118:Yes)、ステップS120へ進む。
ステップS120では、特定機能35Fが、物質のX線透過距離を特定する(ステップS120)。
例えば、特定機能35Fは、記憶回路34に記憶されている、上記ステップS108〜ステップS116の一連の処理毎に算出された誤差の内、最小の誤差を記憶回路34から特定する。そして、特定機能35Fは、特定した最小の誤差に対応する推定距離を、記憶回路34から読取る。これによって、特定機能35Fは、最小の誤差に対応する推定距離を、物質のX線透過距離として特定する。すなわち、特定機能35Fは、弁別対象の物質の各々の、X線透過距離を特定する。
なお、特定機能35Fは、記憶回路34に記憶されている、上記ステップS108〜ステップS116の一連の処理毎に算出された誤差の内、ステップS118の判断において上記一連の処理を予め定めた回数繰り返したと判別したときに算出した誤差を、記憶回路34から特定してもよい、そして、特定機能35Fは、特定した最小の誤差に対応する推定距離を、記憶回路34から読取る。これによって、特定機能35Fは、上記一連の処理を予め定めた回数繰り返したときの推定距離を、物質のX線透過距離として特定してもよい。
また、上記ステップS118の判断において、特定機能35Fが、前回算出した誤差からの変動が予め定めた変動値以下となったか否かを判別することで、ステップS118の判断を行ったと仮定する。この場合、特定機能35Fは、ステップS118で肯定判断したときに算出した誤差に対応する推定距離を、物質のX線透過距離として特定してもよい。このように、特定機能35Fは、算出した誤差の変動が予め定めた変動値以下となる推定距離を、物質のX線透過距離として特定してもよい。
そして、特定機能35Fは、特定したX線透過距離を、対応する物質の種類を示す種類情報に対応付けて、記憶回路34へ記憶する。
次に、生成機能35Gが、ステップS120で記憶回路34へ記憶された、物質の各々のX線透過距離に応じた表示画像を生成する(ステップS122)。次に、表示制御機能35Iが、ステップS122で生成された表示画像を、ディスプレイ32へ表示する(ステップS124)。そして、本ルーチンを終了する。
なお、ステップS122およびステップS124の処理は、ユーザによる入力回路33の操作指示によって表示画像の表示が指示されたときに、実行してもよい。
以上説明したように、本実施の形態のX線CT装置1は、演算機能35Cと、算出機能35Dと、特定機能35Fと、を備える。演算機能35Cは、X線発生器12から被検体Pに照射するX線の照射スペクトル50と、弁別対象の物質のX線透過距離の推定値を示す推定距離と、X線を検出する検出器13の応答特性によって生じるスペクトルの歪みを示す応答情報と、に基づいて、被検体Pを透過したX線の検出器13により検出される推定スペクトル56を演算する。算出機能35Dは、推定スペクトル56と、検出器13で検出された、被検体Pを透過したX線の検出スペクトル54と、の誤差を算出する。特定機能35Fは、誤差に基づいて、物質のX線透過距離を特定する。
このように、本実施の形態のX線CT装置1では、照射スペクトル50から推定距離と応答情報に基づいて演算した推定スペクトル56と、実際に検出されたスペクトルである検出スペクトル54と、の誤差を算出する。そして、X線CT装置1では、誤差に基づいて、物質のX線透過距離を特定する。
このため、本実施の形態のX線CT装置1では、検出スペクトル54から推定スペクトル56を推定する場合などに比べて、推定誤差の影響を抑制することができる。
従って、本実施の形態のX線CT装置1では、高精度の物質弁別を行うことができる。
また、特定機能35Fは、誤差が閾値以下となる推定距離を、物質のX線透過距離として特定することができる。また、特定機能35Fは、誤差が最小となる推定距離を、物質のX線透過距離として特定してもよい。
また、特定機能35Fは、算出した誤差の変動が予め定めた変動値以下となる推定距離を、物質のX線透過距離として特定してもよい。また、特定機能35Fは、推定距離の更新、更新後の推定距離に基づいた推定スペクトル56の演算、および誤差の算出、の一連の処理を予め定めた回数繰り返したときの推定距離を、物質のX線透過距離として特定してもよい。
ここで、演算機能35Cが推定スペクトル56の演算時に最初に用いる推定距離の初期値には、ビームハードニング効果や、検出素子16の応答特性などの影響が含まれている。しかし、本実施の形態のX線CT装置1では、誤差が閾値以下となる推定距離、誤差が最小となる推定距離、誤差の変動が予め定めた変更値以下の推定距離、または、上記一連の処理を予め定めた回数繰り返した時の推定距離を、物質のX線透過距離として特定する。このため、X線CT装置1は、高精度に物質弁別を行うことができる。
また、本実施の形態では、情報処理装置30が用いる検出スペクトル54によって示される、光子カウント数の各々に対応するエネルギー帯のエネルギー幅は、第1のエネルギー幅E1である(図2(D)参照)。また、照射スペクトル50によって示される、光子カウント数の各々に対応するエネルギー帯のエネルギー幅は、第2のエネルギー幅E2である(図2(B)参照)。そして、本実施の形態では、第2のエネルギー幅E2は、第1のエネルギー幅E1以下である。
また、第1のエネルギー幅E1および第2のエネルギー幅E2は、20keV未満や、1keV以下、といった、狭い幅である。
このため、本実施の形態のX線CT装置1では、エネルギー帯内の光子の割合が物質透過中に変化するビームハードニング効果の影響の抑制されたスペクトルを用いて、弁別対象の物質のX線透過距離を特定することができる。
このため、本実施の形態のX線CT装置1では、上記効果に加えて、ビームハードニング効果の影響を低減した、精度の高い物質弁別を行うことができる。
次に、本実施の形態の情報処理装置30のハードウェア構成の一例を説明する。図5は、本実施の形態の情報処理装置30のハードウェア構成図の一例である。情報処理装置30は、CPU(Central Processing Unit)70などの制御装置と、ROM(Read Only Memory)72やRAM(Random Access Memory)74などの記憶装置と、各種機器とのインターフェースであるI/F76と、各部を接続するバス78とを備えており、通常のコンピュータを利用したハードウェア構成となっている。なお、CPU70、ROM72、RAM74、およびI/F76は、バス78を介して互いに接続されている。
上記実施の形態の情報処理装置30では、CPU70が、ROM72からプログラムをRAM74上に読み出して実行することにより、上記各機能がコンピュータ上で実現される。
なお、上記実施の形態の情報処理装置30で実行される上記各処理を実行するためのプログラムは、I/F76を介して接続されたハードディスクに記憶されていてもよい。また、上記実施の形態の情報処理装置30で実行される上記各処理を実行するためのプログラムは、ROM72に予め組み込まれて提供されていてもよい。
また、上記実施の形態の情報処理装置30で実行される上記処理を実行するためのプログラムは、インストール可能な形式または実行可能な形式のファイルでCD−ROM、CD−R、メモリカード、DVD(Digital Versatile Disk)、フレキシブルディスク(FD)等のコンピュータで読み取り可能な記憶媒体に記憶されてコンピュータプログラムプロダクトとして提供されるようにしてもよい。また、上記実施の形態の情報処理装置30で実行される上記処理を実行するためのプログラムを、インターネットなどのネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることにより提供するようにしてもよい。また、上記実施の形態の情報処理装置30で実行される上記処理を実行するためのプログラムを、インターネットなどのネットワーク経由で提供または配布するようにしてもよい。
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。