WO2020250900A1 - 画像処理装置及び画像処理方法、プログラム - Google Patents

Abstract

画像処理装置は、被写体に放射線を照射して撮影を行うことで得られた互いに異なる複数の放射線エネルギーに対応した複数の放射線画像を取得し、取得された複数の放射線画像からエネルギーサブトラクション処理で用いられる複数の減弱率画像を生成する。画像処理装置は、放射線の線量、被写体の厚み、および放射線のエネルギーの少なくともいずれかに依存して生じる減弱率の誤差を低減するように放射線画像または減弱率画像を補正し、補正後の複数の減弱率画像を用いてエネルギーサブトラクション処理により物質特性画像を生成する。

Description

画像処理装置及び画像処理方法、プログラム

　本発明は、画像処理装置及び画像処理方法、プログラムに関する。より具体的には、医療診断における一般撮影などの静止画撮影や透視撮影などの動画撮影に好適に用いられる放射線撮像装置及び放射線撮像システムにおける画像処理装置及び画像処理方法、プログラムに関する。

　現在、Ｘ線による医療画像診断や非破壊検査に用いる撮影装置として、半導体材料によって形成された平面検出器（Ｆｌａｔ　Ｐａｎｅｌ　Ｄｅｔｅｃｔｏｒ、以下ＦＰＤと略す）を用いた放射線撮像装置が普及している。このような放射線撮像装置は、例えば医療画像診断においては、一般撮影のような静止画撮影や、透視撮影のような動画撮影のデジタル撮像装置として用いられている。

　ＦＰＤを用いた撮影方法のひとつに、エネルギーサブトラクションがある。エネルギーサブトラクションでは、まず、管電圧の異なるＸ線を照射するなどして、エネルギーの異なる画像を複数枚取得する。それらを演算することで、骨画像と軟部組織画像に分離するなどの処理を行うことができる（特許文献１）。

特開昭５８－２２１５８０号公報

　エネルギーサブトラクションを用いて骨の厚みと軟部組織の厚みを推定した場合に、厚みの推定値と真値との間に誤差が生じるという課題がある。この結果、骨密度測定の精度が低下したり、推定された骨の厚みや軟部組織の厚みに対して画像処理を適用したときにアーチファクトが生じたりするなどの課題がある。

　本発明は、エネルギーサブトラクションによる物質特性画像の生成において生じる誤差を低減する技術を提供する。

　本発明の一態様による画像処理装置は以下の構成を備える。すなわち、
　被写体に放射線を照射して撮影を行うことで得られた互いに異なる複数の放射線エネルギーに対応した複数の放射線画像を取得する取得手段と、
　前記取得手段により取得された前記複数の放射線画像からエネルギーサブトラクション処理で用いられる複数の減弱率画像を生成する生成手段と、
　前記放射線の線量、前記被写体の厚み、および前記放射線のエネルギーの少なくともいずれかに依存して生じる減弱率の誤差を低減するように、放射線画像または減弱率画像を補正する補正手段と、
　前記補正手段による補正後の前記複数の減弱率画像を用いて前記エネルギーサブトラクション処理により物質特性画像を生成する処理手段と、を備える。

　本願発明によれば、エネルギーサブトラクションによる物質特性画像の生成において生じる誤差が低減される。

　本発明のその他の特徴及び利点は、添付図面を参照とした以下の説明により明らかになるであろう。なお、添付図面においては、同じ若しくは同様の構成には、同じ参照番号を付す。

　添付図面は明細書に含まれ、その一部を構成し、本発明の実施の形態を示し、その記述と共に本発明の原理を説明するために用いられる。
図１は、第１実施形態によるＸ線撮影システムの構成例を示す図である。 図２は、第１実施形態によるＸ線撮像装置の画素等価回路図である。 図３は、第１実施形態によるＸ線撮像装置のタイミングチャートである。 図４は、第１実施形態によるＸ線撮像装置のタイミングチャートである。 図５Ａは、第１実施形態による補正処理を説明する図である。 図５Ｂは、第１実施形態による信号処理のブロック図である。 図５Ｃは、第１実施形態による画像処理のブロック図である。 図６は、減弱率の線量依存を示すグラフ６ａ、画素値の線量依存を示すグラフ６ｂ、第１実施形態による線量依存の補正係数を示すグラフ６ｃを示す図である。 図７Ａは、第１実施形態による補正処理を示す図である。 図７Ｂは、第１実施形態による補正処理を示す図である。 図８は、減弱率の厚み依存を示すグラフ８ａ、第２実施形態に係る厚み依存の補正係数を示すグラフ８ｂ、を示す図である。 図９は、減弱率のスペクトル依存を示すグラフ９ａ、第２実施形態に係るスペクトル依存の補正係数を示すグラフ９ｂ、を示す図である。 図１０は、第２実施形態に係る補正処理を示す図である。 図１１は、第３実施形態によるＸ線撮影のための構成１１ａ、補正値を取得するための構成１１ｂ、および補正値を取得するためのＸ線撮影１１ｃを説明する図である。 図１２は、第３実施形態による物質分離画像の取得処理を説明する図である。 図１３は、第４実施形態によるＸ線撮影のための構成１３ａ、および補正値を取得するためのＸ線撮影の構成１３ｂを示す図である。 図１４は、第５実施形態による物質分離画像の取得処理を示す図である。

　以下、添付図面を参照して実施形態を詳しく説明する。尚、以下の実施形態は特許請求の範囲に係る発明を限定するものではない。実施形態には複数の特徴が記載されているが、これらの複数の特徴の全てが発明に必須のものとは限らず、また、複数の特徴は任意に組み合わせられてもよい。さらに、添付図面においては、同一若しくは同様の構成に同一の参照番号を付し、重複した説明は省略する。

　なお、本発明における放射線には、放射線崩壊によって放出される粒子（光子を含む）の作るビームであるα線、β線、γ線などの他に、同程度以上のエネルギーを有するビーム、例えばＸ線や粒子線、宇宙線なども、含まれるものとする。以下の実施形態では、放射線の一例としてＸ線を用いた装置を説明する。したがって、以下では、放射線、放射線画像、放射線エネルギー、放射線スペクトル、放射線量、放射線発生装置、放射線撮像装置、放射線撮像システムとして、それぞれＸ線、Ｘ線画像、Ｘ線エネルギー、Ｘ線スペクトル、Ｘ線量、Ｘ線発生装置、Ｘ線撮像装置、Ｘ線撮像システムとして説明する。

　＜第１実施形態＞
　図１は、第１実施形態に係る、放射線撮像システムの一例としてのＸ線撮像システムの構成例を示すブロック図である。第１実施形態のＸ線撮像システムは、Ｘ線発生装置１０１、Ｘ線制御装置１０２、撮像制御装置１０３、Ｘ線撮像装置１０４を備える。

　Ｘ線発生装置１０１は、Ｘ線を発生し、被写体を曝射する。Ｘ線制御装置１０２は、Ｘ線発生装置１０１におけるＸ線の発生を制御する。撮像制御装置１０３は、例えば、１つまたは複数のプロセッサー（ＣＰＵ）とメモリを有し、プロセッサーがメモリに格納されたプログラムを実行してＸ線画像の取得及び画像処理を行う。なお、撮像制御装置１０３による画像処理を含む各処理は、専用のハードウエアにより実現されてもよいし、ハードウエアとソフトウエアの協働により実現されてもよい。Ｘ線撮像装置１０４は、Ｘ線を可視光に変換する蛍光体１０５と、可視光を検出する二次元検出器１０６を有する。二次元検出器は、Ｘ線量子を検出する画素２０をＸ列×Ｙ行のアレイ状に配置したセンサであり、画像情報を出力する。

　撮像制御装置１０３は、上述したプロセッサーにより放射線画像を処理する画像処理装置として機能する。取得部１３１、補正部１３２、信号処理部１３３、画像処理部１３４は、画像処理装置としての機能構成例を示している。取得部１３１は、被写体に放射線を照射して撮影を行うことで得られた互いに異なる複数の放射線エネルギーに対応した複数の放射線画像を取得する。補正部１３２は、取得部１３１により取得された複数の放射線画像を補正してエネルギーサブトラクション処理で用いられる複数の画像を生成する。また、補正部１３２は、線量に依存して生じる誤差、被写体の厚みに依存して生じる誤差、および放射線のエネルギーに依存して生じる誤差のうちの少なくとも何れかを低減するための補正を行う。補正部１３２の詳細については後述する。信号処理部１３３は、補正部１３２により生成された複数の画像を用いて物質特性画像を生成する。物質特性画像とは、例えば骨と軟部組織というように物質を分離して表す物質分離画像、実効原子番号とその面密度を表す物質識別画像など、エネルギーサブトラクション処理において取得される画像である。信号処理部１３３の詳細は後述する。画像処理部１３４は、取得された物質特性画像を用いて、仮想単色Ｘ線画像を生成する。画像処理部１３４の詳細は後述する。

　図２は、第１実施形態に係る画素２０の等価回路図である。画素２０は、光電変換素子２０１と、出力回路部２０２とを含む。光電変換素子２０１は、典型的にはフォトダイオードでありうる。出力回路部２０２は、増幅回路部２０４、クランプ回路部２０６、サンプルホールド回路部２０７、選択回路部２０８を含む。

　光電変換素子２０１は、電荷蓄積部を含み、該電荷蓄積部は、増幅回路部２０４のＭＯＳトランジスタ２０４ａのゲートに接続されている。ＭＯＳトランジスタ２０４ａのソースは、ＭＯＳトランジスタ２０４ｂを介して電流源２０４ｃに接続されている。ＭＯＳトランジスタ２０４ａと電流源２０４ｃとによってソースフォロア回路が構成されている。ＭＯＳトランジスタ２０４ｂは、そのゲートに供給されるイネーブル信号ＥＮがアクティブレベルになるとオンしてソースフォロア回路を動作状態にするイネーブルスイッチである。

　図２に示す例では、光電変換素子２０１の電荷蓄積部およびＭＯＳトランジスタ２０４ａのゲートが共通のノードを構成していて、このノードは、電荷蓄積部に蓄積された電荷を電圧に変換する電荷電圧変換部として機能する。即ち、電荷電圧変換部には、電荷蓄積部に蓄積された電荷Ｑと電荷電圧変換部が有する容量値Ｃとによって定まる電圧Ｖ（＝Ｑ／Ｃ）が現れる。電荷電圧変換部は、リセットスイッチ２０３を介してリセット電位Ｖｒｅｓに接続されている。リセット信号ＰＲＥＳがアクティブレベルになると、リセットスイッチ２０３がオンして、電荷電圧変換部の電位がリセット電位Ｖｒｅｓにリセットされる。

　クランプ回路部２０６は、リセットした電荷電圧変換部の電位に応じて増幅回路部２０４によって出力されるノイズをクランプ容量２０６ａによってクランプする。つまり、クランプ回路部２０６は、光電変換素子２０１で光電変換により発生した電荷に応じてソースフォロア回路から出力された信号から、このノイズをキャンセルするための回路である。このノイズはリセット時のｋＴＣノイズを含む。クランプは、クランプ信号ＰＣＬをアクティブレベルにしてＭＯＳトランジスタ２０６ｂをオン状態にした後に、クランプ信号ＰＣＬを非アクティブレベルにしてＭＯＳトランジスタ２０６ｂをオフ状態にすることによってなされる。クランプ容量２０６ａの出力側は、ＭＯＳトランジスタ２０６ｃのゲートに接続されている。ＭＯＳトランジスタ２０６ｃのソースは、ＭＯＳトランジスタ２０６ｄを介して電流源２０６ｅに接続されている。ＭＯＳトランジスタ２０６ｃと電流源２０６ｅとによってソースフォロア回路が構成されている。ＭＯＳトランジスタ２０６ｄは、そのゲートに供給されるイネーブル信号ＥＮ０がアクティブレベルになるとオンしてソースフォロア回路を動作状態にするイネーブルスイッチである。

　光電変換素子２０１で光電変換により発生した電荷に応じてクランプ回路部２０６から出力される信号は、光信号として、光信号サンプリング信号ＴＳがアクティブレベルになることによってスイッチ２０７Ｓａを介して容量２０７Ｓｂに書き込まれる。電荷電圧変換部の電位をリセットした直後にＭＯＳトランジスタ２０６ｂをオン状態とした際にクランプ回路部２０６から出力される信号は、クランプ電圧である。ノイズ信号は、ノイズサンプリング信号ＴＮがアクティブレベルになることによってスイッチ２０７Ｎａを介して容量２０７Ｎｂに書き込まれる。このノイズ信号には、クランプ回路部２０６のオフセット成分が含まれる。スイッチ２０７Ｓａと容量２０７Ｓｂによって信号サンプルホールド回路２０７Ｓが構成され、スイッチ２０７Ｎａと容量２０７Ｎｂによってノイズサンプルホールド回路２０７Ｎが構成される。サンプルホールド回路部２０７は、信号サンプルホールド回路２０７Ｓとノイズサンプルホールド回路２０７Ｎとを含む。

　駆動回路部が行選択信号をアクティブレベルに駆動すると、容量２０７Ｓｂに保持された信号（光信号）がＭＯＳトランジスタ２０８Ｓａおよび行選択スイッチ２０８Ｓｂを介して信号線２１Ｓに出力される。また、同時に、容量２０７Ｎｂに保持された信号（ノイズ）がＭＯＳトランジスタ２０８Ｎａおよび行選択スイッチ２０８Ｎｂを介して信号線２１Ｎに出力される。ＭＯＳトランジスタ２０８Ｓａは、信号線２１Ｓに設けられた不図示の定電流源とソースフォロア回路を構成する。同様に、ＭＯＳトランジスタ２０８Ｎａは、信号線２１Ｎに設けられた不図示の定電流源とソースフォロア回路を構成する。ＭＯＳトランジスタ２０８Ｓａと行選択スイッチ２０８Ｓｂによって信号用選択回路部２０８Ｓが構成され、ＭＯＳトランジスタ２０８Ｎａと行選択スイッチ２０８Ｎｂによってノイズ用選択回路部２０８Ｎが構成される。選択回路部２０８は、信号用選択回路部２０８Ｓとノイズ用選択回路部２０８Ｎとを含む。

　画素２０は、隣接する複数の画素２０の光信号を加算する加算スイッチ２０９Ｓを有してもよい。加算モード時には、加算モード信号ＡＤＤがアクティブレベルになり、加算スイッチ２０９Ｓがオン状態になる。これにより、隣接する画素２０の容量２０７Ｓｂが加算スイッチ２０９Ｓによって相互に接続されて、光信号が平均化される。同様に、画素２０は、隣接する複数の画素２０のノイズを加算する加算スイッチ２０９Ｎを有してもよい。加算スイッチ２０９Ｎがオン状態になると、隣接する画素２０の容量２０７Ｎｂが加算スイッチ２０９Ｎによって相互に接続されて、ノイズが平均化される。加算部２０９は、加算スイッチ２０９Ｓと加算スイッチ２０９Ｎを含む。

　また、画素２０は、感度を変更するための感度変更部２０５を有してもよい。画素２０は、例えば、第１感度変更スイッチ２０５ａおよび第２感度変更スイッチ２０５'ａ、並びにそれらに付随する回路素子を含みうる。第１変更信号ＷＩＤＥがアクティブレベルになると、第１感度変更スイッチ２０５ａがオンして、電荷電圧変換部の容量値に第１付加容量２０５ｂの容量値が追加される。これによって画素２０の感度が低下する。第２変更信号ＷＩＤＥ２がアクティブレベルになると、第２感度変更スイッチ２０５'ａがオンして、電荷電圧変換部の容量値に第２付加容量２０５'ｂの容量値が追加される。これによって画素２０１の感度が更に低下する。このように画素２０の感度を低下させる機能を追加することによって、より大きな光量を受光することが可能となり、ダイナミックレンジを広げることができる。第１変更信号ＷＩＤＥがアクティブレベルになる場合には、イネーブル信号ＥＮｗをアクティブレベルにして、ＭＯＳトランジスタ２０４ａに変えてＭＯＳトランジスタ２０４'ａをソースフォロア動作させてもよい。

　Ｘ線撮像装置１０４は、二次元検出器１０６から以上のような画素回路の出力を読み出し、不図示のＡＤ変換器でデジタル値に変換した後、撮像制御装置１０３に画像を転送する。

　次に、上述した構成を備えた第１実施形態のＸ線撮像システムの動作について説明する。図３は、第１実施形態に係るＸ線撮像システムにおいてエネルギーサブトラクションに提供するための、互いにエネルギーの異なる複数のＸ線画像を得る場合のＸ線撮像装置１０４の駆動タイミングを示す。図３中の波形は横軸を時間として、Ｘ線の曝射、同期信号、光電変換素子２０１のリセット、サンプルホールド回路２０７、信号線２１からの画像の読み出しのタイミングを示している。

　リセット信号により光電変換素子２０１のリセットが行われてからＸ線が曝射される。Ｘ線の管電圧は理想的には矩形波となるが、管電圧の立ち上がりと立下りには有限の時間がかかる。特に、パルスＸ線で曝射時間が短い場合は、管電圧はもはや矩形波とはみなせず、Ｘ線３０１～３０３に示すような波形となる。立ち上がり期のＸ線３０１、安定期のＸ線３０２、立下り期のＸ線３０３ではそれぞれＸ線のエネルギーが異なる。したがって、サンプルホールドによって区切られる期間の放射線に対応したＸ線画像を得ることにより、互いにエネルギーが異なる複数種類のＸ線画像が得られる。

　Ｘ線撮像装置１０４は、立ち上がり期のＸ線３０１が曝射された後に、ノイズサンプルホールド回路２０７Ｎでサンプリングを行い、さらに安定期のＸ線３０２が曝射された後に信号サンプルホールド回路２０７Ｓでサンプリングを行う。その後、Ｘ線撮像装置１０４は、信号線２１Ｎと信号線２１Ｓの差分を画像として読み出す。このとき、ノイズサンプルホールド回路２０７Ｎには立ち上がり期のＸ線３０１の信号（Ｒ_１）が保持され、信号サンプルホールド回路２０７Ｓには立ち上がり期のＸ線３０１の信号と安定期のＸ線３０２の信号（Ｂ）の和（Ｒ_１＋Ｂ）が保持されている。従って、安定期のＸ線３０２の信号に対応した画像３０４が読み出される。

　次に、Ｘ線撮像装置１０４は、立下り期のＸ線３０３の曝射と、画像３０４の読み出しとが完了してから、再び信号サンプルホールド回路２０７Ｓでサンプリングを行う。その後、Ｘ線撮像装置１０４は、光電変換素子２０１のリセットを行い、再びノイズサンプルホールド回路２０７Ｎでサンプリングを行い、信号線２１Ｎと信号線２１Ｓの差分を画像として読み出す。このとき、ノイズサンプルホールド回路２０７ＮにはＸ線が曝射されていない状態の信号が保持され、信号サンプルホールド回路２０７Ｓには立ち上がり期のＸ線３０１の信号と安定期のＸ線３０２と立下り期のＸ線３０３の信号（Ｒ_２）の和（Ｒ_１＋Ｂ＋Ｒ_２）が保持されている。従って、立ち上がり期のＸ線３０１の信号と安定期のＸ線３０２の信号と立下り期のＸ線３０３の信号に対応した画像３０６が読み出される。その後、画像３０６と画像３０４の差分を計算することで、立ち上がり期のＸ線３０１と立下り期のＸ線３０３の和に対応した画像３０５が得られる。この計算は、Ｘ線撮像装置１０４で行われてもよいし、撮像制御装置１０３で行われてもよい。

　サンプルホールド回路２０７及び光電変換素子２０１のリセットを行うタイミングは、Ｘ線発生装置１０１からＸ線の曝射が開始されたことを示す同期信号３０７を用いて決定される。Ｘ線の曝射開始を検出する方法としては、Ｘ線発生装置１０１の管電流を測定し、電流値が予め設定された閾値を上回るか否かを判定する構成を用いることができるがこれに限られるものではない。例えば、光電変換素子２０１のリセットが完了した後、画素２０を繰り返して読み出し、画素値が予め設定された閾値を上回るか否かを判定することによりＸ線の曝射開始を検出する構成が用いられてもよい。あるいは、例えば、Ｘ線撮像装置１０４に二次元検出器１０６とは異なるＸ線検出器を内蔵し、その測定値が予め設定された閾値を上回るか否かを判定することによりＸ線の曝射開始を検出する構成が用いられてもよい。いずれの方式の場合も、Ｘ線の曝射開始を示す同期信号３０７の入力から予め指定した時間が経過した後に、信号サンプルホールド回路２０７Ｓのサンプリング、ノイズサンプルホールド回路２０７Ｎのサンプリング、光電変換素子２０１のリセットが行われる。

　以上のようにして、パルスＸ線の安定期に対応した画像３０４と、立ち上がり期と立下り期の和に対応した画像３０５が得られる。これら二枚のＸ線画像を形成する際に曝射されたＸ線のエネルギーは互いに異なるため、これらＸ線画像間で演算を行うことでエネルギーサブトラクション処理を行うことができる。

　図４は、第１実施形態に係るＸ線撮像システムにおいてエネルギーサブトラクションに提供するための、互いにエネルギーの異なる複数のＸ線画像を得る、図３とは異なるＸ線撮像装置１０４の駆動タイミングを示す。図３とは、Ｘ線発生装置１０１の管電圧を能動的に切り替えている点で異なる。

　まず、光電変換素子２０１のリセットが行われた後、Ｘ線発生装置１０１は低エネルギーのＸ線４０１の曝射を行う。この状態で、Ｘ線撮像装置１０４は、ノイズサンプルホールド回路２０７Ｎによりサンプリングを行う。その後、Ｘ線発生装置１０１は、管電圧を切り替えて高エネルギーのＸ線４０２の曝射を行う。この状態で、Ｘ線撮像装置１０４は、信号サンプルホールド回路２０７Ｓによりサンプリングを行う。その後、Ｘ線発生装置１０１は、管電圧を切り替えて低エネルギーのＸ線４０３の曝射を行う。Ｘ線撮像装置１０４は、信号線２１Ｎと信号線２１Ｓの差分を画像として読み出す。このとき、ノイズサンプルホールド回路２０７Ｎには低エネルギーのＸ線４０１の信号（Ｒ_１）が保持され、信号サンプルホールド回路２０７Ｓには低エネルギーのＸ線４０１の信号と高エネルギーのＸ線４０２の信号（Ｂ）の和（Ｒ_１＋Ｂ）が保持されている。従って、高エネルギーのＸ線４０２の信号に対応した画像４０４が読み出される。

　次に、Ｘ線撮像装置１０４は、低エネルギーのＸ線４０３の曝射と、画像４０４の読み出しとが完了してから、再び信号サンプルホールド回路２０７Ｓでサンプリングを行う。その後、Ｘ線撮像装置１０４は、光電変換素子２０１のリセットを行い、再びノイズサンプルホールド回路２０７Ｎでサンプリングを行い、信号線２１Ｎと信号線２１Ｓの差分を画像として読み出す。このとき、ノイズサンプルホールド回路２０７ＮにはＸ線が曝射されていない状態の信号が保持され、信号サンプルホールド回路２０７Ｓには低エネルギーのＸ線４０１の信号と高エネルギーのＸ線４０２と低エネルギーのＸ線４０３の信号（Ｒ_２）の和（Ｒ_１＋Ｂ＋Ｒ_２）が保持されている。従って、低エネルギーのＸ線４０１の信号と高エネルギーのＸ線４０２の信号と低エネルギーのＸ線４０３の信号に対応した画像４０６が読み出される。

　その後、画像４０６と画像４０４の差分を計算することで、低エネルギーのＸ線４０１と低エネルギーのＸ線４０３の和に対応した画像４０５が得られる。この計算は、Ｘ線撮像装置１０４で行われてもよいし、撮像制御装置１０３で行われてもよい。同期信号４０７については、図３と同様である。このように、管電圧を能動的に切り替えながら画像を取得することで、図３の方法に比べて低エネルギーと高エネルギーの放射線画像の間のエネルギー差をより大きくすることが出来る。

　次に、撮像制御装置１０３によるエネルギーサブトラクション処理について説明する。第１実施形態におけるエネルギーサブトラクション処理は、補正部１３２による補正処理、信号処理部１３３による信号処理、画像処理部１３４による画像処理の３段階に分かれている。以下、それぞれの処理について説明する。

　（補正処理の説明）
　補正処理は、Ｘ線撮像装置１０４から取得された複数の放射線画像を処理してエネルギーサブトラクション処理における後述の信号処理で用いられる複数の画像を生成する処理である。図５Ａに、第１実施形態に係るエネルギーサブトラクション処理のための補正処理を示す。まず、取得部１３１は、Ｘ線撮像装置１０４にＸ線を曝射しない状態での撮像を行わせ、図３または図４に示した駆動で画像を取得する。この駆動により２枚の画像が読み出される。以下、１枚目の画像（画像３０４または画像４０４）をＦ＿ＯＤＤ、２枚目の画像（画像３０６または画像４０６）をＦ＿ＥＶＥＮとする。Ｆ＿ＯＤＤとＦ＿ＥＶＥＮは、Ｘ線撮像装置１０４の固定パターンノイズ（ＦＰＮ：Ｆｉｘｅｄ　Ｐａｔｔｅｒｎ　Ｎｏｉｓｅ）に対応する画像である。

　次に、取得部１３１は、被写体がない状態でＸ線撮像装置１０４にＸ線を曝射して撮像を行わせ、図３又は図４に示した駆動によりＸ線撮像装置１０４から出力されるゲイン補正用の画像を取得する。この駆動により、上記と同様に２枚の画像が読み出される。以下、１枚目のゲイン補正用の画像（画像３０４または画像４０４）をＷ＿ＯＤＤ、２枚目のゲイン補正用の画像（画像３０６または画像４０６）をＷ＿ＥＶＥＮとする。Ｗ＿ＯＤＤとＷ＿ＥＶＥＮは、Ｘ線撮像装置１０４のＦＰＮとＸ線による信号の和に対応する画像である。補正部１３２は、Ｗ＿ＯＤＤからＦ＿ＯＤＤを、Ｗ＿ＥＶＥＮからＦ＿ＥＶＥＮを減算することで、Ｘ線撮像装置１０４のＦＰＮが除去された画像ＷＦ＿ＯＤＤとＷＦ＿ＥＶＥＮを得る。これをオフセット補正と呼ぶ。

　ＷＦ＿ＯＤＤは安定期のＸ線３０２に対応する画像であり、ＷＦ＿ＥＶＥＮは立ち上がり期のＸ線３０１、安定期のＸ線３０２、立下り期のＸ線３０３の和に対応する画像である。従って、補正部１３２は、ＷＦ＿ＥＶＥＮからＷＦ＿ＯＤＤを減算することで、立ち上がり期のＸ線３０１と立下り期のＸ線３０３の和に対応する画像を得る。このように、複数枚の画像の減算により、サンプルホールドによって区切られる特定の期間のＸ線に対応した画像を得る処理を色補正と呼ぶ。立ち上がり期のＸ線３０１と立下り期のＸ線３０３のエネルギーは、安定期のＸ線３０２のエネルギーに比べて低い。従って、色補正により、ＷＦ＿ＥＶＥＮからＷＦ＿ＯＤＤを減算することで、被写体がない場合の低エネルギー画像Ｗ＿Ｌｏｗが得られる。また、ＷＦ＿ＯＤＤから、被写体がない場合の高エネルギー画像Ｗ＿Ｈｉｇｈが得られる。

　次に、取得部１３１は、被写体がある状態でＸ線撮像装置１０４にＸ線を曝射して撮像を行わせ、図３または図４に示した駆動によりＸ線撮像装置１０４から出力される画像を取得する。このとき２枚の画像が読み出される。以下、１枚目の画像（画像３０４または画像４０４）をＸ＿ＯＤＤ、２枚目の画像（画像３０６または画像４０６）をＸ＿ＥＶＥＮとする。補正部１３２は、被写体がない場合と同様のオフセット補正および色補正を行うことで、被写体がある場合の低エネルギー画像Ｘ＿Ｌｏｗと、被写体がある場合の高エネルギー画像Ｘ＿Ｈｉｇｈを得る。

　ここで、被写体の厚みをｄ、被写体の線減弱係数をμ、被写体がない場合の画素２０の出力をＩ_０、被写体がある場合の画素２０の出力をＩとすると、以下の［数１］式が成り立つ。

　［数１］を変形すると、以下の［数２］式が得られる。［数２］の右辺は被写体の減弱率を示す。被写体の減弱率は０～１の間の実数である。

　従って、補正部１３２は、被写体がある場合の低エネルギー画像Ｘ＿Ｌｏｗを、被写体がない場合の低エネルギー画像Ｗ＿Ｌｏｗで除算することで、低エネルギーにおける減弱率の画像Ｌを得る。同様に、補正部１３２は、被写体がある場合の高エネルギー画像Ｘ＿Ｈｉｇｈを、被写体がない場合の高エネルギー画像Ｗ＿Ｈｉｇｈで除算することで、高エネルギーにおける減弱率の画像Ｈを得る。このように、被写体ありの状態で得られた放射線画像に基づいて得られた画像を被写体なしの状態で得られた放射線画像に基づいて得られた画像で除算することにより減弱率の画像を取得する処理をゲイン補正と呼ぶ。以上が、第１実施形態の補正部１３２による補正処理（但し、後述の、線量依存補正を除く）の説明である。

　（信号処理の説明）
　図５Ｂに、第１実施形態に係るエネルギーサブトラクション処理の信号処理のブロック図を示す。信号処理部１３３は、補正部１３２から得られる複数の画像を用いて物質特性画像を生成する。以下では、骨の厚さの画像Ｂと軟部組織の厚さの画像Ｓからなる物質分離画像の生成を説明する。信号処理部１３３は、以下の処理により図５Ａに示した補正によって得られた低エネルギーにおける減弱率の画像Ｌと高エネルギーにおける減弱率の画像Ｈから、骨の厚さの画像Ｂと軟部組織の厚さの画像Ｓを求める。

　まず、Ｘ線フォトンのエネルギーをＥ、エネルギーＥにおけるフォトン数をＮ（Ｅ）、骨の厚さをＢ、軟部組織の厚さをＳ、エネルギーＥにおける骨の線減弱係数をμ_Ｂ（Ｅ）、エネルギーＥにおける軟部組織の線減弱係数をμ_Ｓ（Ｅ）、減弱率をＩ／Ｉ_０とすると、以下の［数３］式が成り立つ。

　エネルギーＥにおけるフォトン数Ｎ（Ｅ）は、Ｘ線のスペクトルである。Ｘ線のスペクトルは、シミュレーション又は実測により得られる。また、エネルギーＥおける骨の線減弱係数μ_Ｂ（Ｅ）とエネルギーＥおける軟部組織の線減弱係数μ_Ｓ（Ｅ）は、それぞれＮＩＳＴ（Ｎａｔｉｏｎａｌ　Ｉｎｓｔｉｔｕｔｅ　ｏｆ　Ｓｔａｎｄａｒｄｓ　ａｎｄ　Ｔｅｃｈｎｏｌｏｇｙ）などのデータベースから得られる。したがって、［数３］によれば、任意の骨の厚さＢ、軟部組織の厚さＳ、Ｘ線のスペクトルＮ（Ｅ）における減弱率Ｉ／Ｉ_０を計算することが可能である。

　ここで、低エネルギーのＸ線におけるスペクトルをＮ_Ｌ（Ｅ）、高エネルギーのＸ線におけるスペクトルをＮ_Ｈ（Ｅ）とすると、以下の［数４］の各式が成り立つ。なお、Ｌは低エネルギーの減弱率の画像における画素値、Ｈは高エネルギーの減弱率の画像における画素値である。

　［数４］の非線形連立方程式を解くことで、骨の厚みＢと軟部組織の厚みＳが求まる。非線形連立方程式を解く代表的な方法として、ここではニュートンラフソン法を用いた場合について説明する。まず、ニュートンラフソン法の反復回数をｍ、ｍ回目の反復後の骨の厚みをＢ^ｍ、ｍ回目の反復後の軟部組織の厚みをＳ^ｍとしたとき、ｍ回目の反復後の高エネルギーの減弱率Ｈ^ｍ、ｍ回目の反復後の低エネルギーの減弱率Ｌ^ｍは以下の［数５］で表される。

　また、厚みが微小に変化したときの減弱率の変化率を、以下の［数６］で表す。

　このとき、ｍ＋１回目の反復後の骨の厚みＢ^ｍ＋１と軟部組織の厚みＳ^ｍ＋１を、高エネルギーの減弱率Ｈと低エネルギーの減弱率Ｌを用いて、以下の［数７］で表す。

　２ｘ２の行列の逆行列は、行列式をｄｅｔとすると、クラメルの公式より以下の［数８］で表される。

　従って、［数７］に［数８］を代入すると、以下の[数９]が求まる。

　以上のような計算を繰り返すことで、ｍ回目の反復後の高エネルギーの減弱率Ｈ^ｍと実測した高エネルギーの減弱率Ｈの差分が限りなく０に近づいていく。低エネルギーの減弱率Ｌについても同様である。これによって、ｍ回目の反復後の骨の厚みＢ^ｍが骨の厚みＢに収束し、ｍ回目の軟部組織の厚みＳ^ｍが軟部組織の厚みＳに収束する。以上のようにして、［数４］に示した非線形連立方程式を解くことができる。従って、全ての画素について［数４］を計算することで、低エネルギーにおける減弱率の画像Ｌと高エネルギーにおける減弱率の画像Ｈから、骨の厚さの画像Ｂ、軟部組織の厚さの画像Ｓを得ることができる。

　なお、第１実施形態では、骨の厚さＢと軟部組織の厚さＳを算出していたが、本発明はこのような形態に限定されない。例えば、水の厚さＷと造影剤の厚さＩを算出してもよい。すなわち、任意の二種類の物質の厚さに分解してもよい。また、図５Ａに示した補正によって得られた低エネルギーにおける減弱率の画像Ｌと高エネルギーにおける減弱率の画像Ｈから、実効原子番号Ｚの画像と面密度Ｄの画像を求めてもよい。実効原子番号Ｚとは混合物の等価的な原子番号のことであり、面密度Ｄとは被写体の密度［ｇ／ｃｍ^３］と被写体の厚み［ｃｍ］の積である。

　また、第１実施形態では、ニュートンラフソン法を用いて非線形連立方程式を解いていた。しかしながら本発明はこのような形態に限定されない。例えば、最小二乗法や二分法などの反復解法を用いてもよい。また、第１実施形態では非線形連立方程式を反復解法で解いていたが、本発明はこのような形態に限定されない。様々な組み合わせの高エネルギーの減弱率Ｈと低エネルギーの減弱率Ｌに対する骨の厚みＢや軟部組織の厚みＳを事前に求めてテーブルを生成し、このテーブルを参照することで骨の厚みＢや軟部組織の厚みＳを高速に求める構成を用いても良い。

　（画像処理の説明）
　図５Ｃに、第１実施形態に係るエネルギーサブトラクション処理の画像処理のブロック図を示す。第１実施形態の画像処理部１３４は、図５Ｂに示した信号処理によって得られた骨の厚さの画像Ｂ、軟部組織の厚さの画像Ｓから、仮想単色Ｘ線画像を得る画像処理を行う。仮想単色Ｘ線画像とは、単一のエネルギーのＸ線を照射した場合に得られることが想定される画像のことである。例えば、仮想単色Ｘ線のエネルギーをＥ_Ｖとしたとき、仮想単色Ｘ線画像Ｖは以下の［数１０］で求められる。

　仮想単色Ｘ線画像は、エネルギーサブトラクションと三次元再構成を組み合わせたＤｕａｌ　Ｅｎｅｒｇｙ　ＣＴで利用されている。このとき、仮想単色Ｘ線画像のＣＮＲ（Ｃｏｎｔｒａｓｔ－ｔｏ－Ｎｏｉｓｅ　Ｒａｔｉｏ）を向上させるために、仮想単色Ｘ線のエネルギーＥ_Ｖを変更する。例えば、骨の線減弱係数μ_Ｂ（Ｅ）は、軟部組織の線減弱係数μ_Ｓ（Ｅ）に比べて大きい。しかしながら、仮想単色Ｘ線のエネルギーＥ_Ｖが大きくなるほどその差は小さくなる。従って、骨画像のノイズによる仮想単色Ｘ線画像のノイズ増加が抑制される。一方で、仮想単色Ｘ線のエネルギーＥ_Ｖが小さくなるほど、μ_Ｂ（Ｅ）とμ_Ｓ（Ｅ）の差が大きくなるため、仮想単色Ｘ線画像のコントラストが大きくなる。すなわち、仮想単色Ｘ線画像のエネルギーＥ_Ｖには適切な値が存在する。

　なお、第１実施形態では骨の厚さＢと軟部組織の厚さＳから仮想単色Ｘ線画像を生成していたが、本発明はこのような形態に限定されない。上述したように、実効原子番号Ｚと面密度Ｄを算出してから、実効原子番号Ｚと面密度Ｄを用いて仮想単色Ｘ線画像を生成してもよい。また、複数のエネルギーＥ_Ｖで生成した複数の仮想単色Ｘ線画像を合成することで、合成Ｘ線画像を生成してもよい。合成Ｘ線画像とは、任意のスペクトルのＸ線を照射した場合に得られることが想定される画像のことである。

　また、第１実施形態の画像処理では仮想単色Ｘ線画像を生成していたが、本発明はこのような形態に限定されない。骨の厚さの画像Ｂや軟部組織の厚さの画像Ｓをそのまま表示してもよい。また、骨の厚さの画像Ｂや軟部組織の厚さの画像Ｓに、リカーシブフィルタ等の時間方向のフィルタや、ガウシアンフィルタ等の空間方向のフィルタをかけるなどして得られた画像を表示するようにしてもよい。また、造影剤を注入する前後の低エネルギーの画像（減弱率）と高エネルギーの画像（減弱率）を用いて、骨のＤＳＡ画像（Ｄｉｇｉｔａｌ　Ｓｕｂｔｒａｃｔｉｏｎ　Ａｎｇｉｏｇｒａｐｈｙ）を得て、これを表示するようにしてもよい。すなわち、本実施形態における画像処理とは、信号処理後の画像に対して任意の演算を行う処理であると言える。

　なお、ＤＳＡ画像は、例えば次のようにして取得される。まず、造影剤を注入する前に、Ｘ線撮影を行って低エネルギーにおける減弱率の画像Ｌ_Ｍと高エネルギーにおける減弱率の画像Ｈ_Ｍを得る。そして、画像Ｌ_Ｍと画像Ｈ_Ｍから骨の厚さのマスク画像Ｂ_Ｍと軟部組織の厚さのマスク画像Ｓ_Ｍを求める。次に、造影剤を注入した後に撮影した低エネルギーにおける減弱率のライブ画像Ｌ_Ｌと高エネルギーにおける減弱率のライブ画像Ｈ_Ｌから、骨の厚さのライブ画像Ｂ_Ｌと軟部組織の厚さのライブ画像Ｓ_Ｌを求める。骨の厚さのライブ画像Ｂ_Ｌから骨の厚さのマスク画像Ｂ_Ｍを引くことで骨のＤＳＡ画像Ｂ_ＤＳＡが得られる。

　以上、本実施形態の補正処理、信号処理、画像処理について説明した。次に、補正部１３２により実行される、推定値に生じる誤差を低減する処理について説明する。図５Ａ～図５Ｃを用いて説明したように、第１実施形態におけるエネルギーサブトラクション処理は、補正処理、信号処理、画像処理の３ステップで構成される。このとき、［数４］を解くことで求めた骨の厚みＢや軟部組織の厚みＳを、厚みの推定値とする。また、定規等で計測した厚みを、厚みの真値とする。補正と信号処理が適切に行われれば、厚みの推定値と厚みの真値は一致するはずである。しかしながら本発明者らが検討を行ったところ、上記のエネルギーサブトラクション処理により得られた厚みの推定値は、厚みの真値に必ずしも一致しないことが判明した。厚みの推定値と真値の誤差が大きくなると、画像処理後の画像にアーチファクトが生じるなどの課題がある。

　本発明者らが検討を行った結果、誤差の原因として、散乱線、減弱率の線量依存、減弱率の厚み依存と減弱率のスペクトル依存などがあることが判明した。本実施形態においては、説明を簡略化するため、散乱線については除去または補正済であるものとし、それ以外の原因による誤差を低減する方法を説明する。すなわち、補正部１３２は、線量に依存して生じる誤差、被写体の厚みに依存して生じる誤差、および放射線のエネルギーに依存して生じる誤差のうちの少なくとも何れかを低減するように、複数の画像を生成する過程（図５Ａ）で得られる画像の各画素の画素値または減弱率を補正する処理を含む。なお、第１実施形態では減弱率の線量依存に係る補正を説明する。減弱率の厚みおよびスペクトルへの依存に係る補正については、第２実施形態で説明する。

　図６のグラフ６ａに、第１実施形態に係る減弱率の線量依存を示す。減弱率Ｉ／Ｉ_０は、［数３］で表される。ここで、Ｘ線発生装置１０１の管電圧を一定にしたまま管電流や曝射時間を調節する、あるいはＸ線発生装置１０１とＸ線撮像装置１０４の距離を調節するなどして、スペクトルＮ（Ｅ）をα倍（０＜α）にする。この場合、［数３］の右辺の分母と分子が共にα倍になって相殺される。すなわち、減弱率Ｉ／Ｉ_０は、線量が変化しても一定になると考えられる。しかしながら、管電圧やフィルタ等を変えず、すなわちスペクトルを変えずに線量のみを変化させながら減弱率を測定したところ、グラフ６ａのように、減弱率が線量の変化に対して一定にならないことが判明した。

　図６のグラフ６ｂに、第１実施形態に係るＸ線撮像装置１０４から得られる画素値の線量依存を示す。Ｘ線撮像装置１０４の画素２０から読み出された画素値は、理想的には線量に比例した値となるはずである。しかしながら、管電圧やフィルタ等を変えず、すなわちスペクトルを変えずに線量を変化させながら画素値を測定したところ、グラフ６ｂのように、画素値の測定値が理想の直線から外れることが判明した。このようなＸ線撮像装置１０４の積分非直線性誤差（ＩＮＬ：Ｉｎｔｅｇｒａｌ　Ｎｏｎ－Ｌｉｎｅａｒｉｔｙ）は、グラフ６ａのように減弱率が一定にならないことの一因である。

　第１実施形態では、補正部１３２が、グラフ６ｂにおける測定値が理想の直線（所定の直線）に近づくように補正を行うことにより、線量の変化に対して減弱率が一定になるようにする。以下、このような補正を線量依存補正と呼ぶ。図６のグラフ６ｃに、第１実施形態に係る線量依存補正の補正係数を示す。第１実施形態では、画素の測定値に対してグラフ６ｃのような補正係数をかけることで、線量依存補正によって、グラフ６ｂの測定値が理想の直線に近づき、グラフ６ａにおける減弱率が一定に近づく。例えば、画素の測定値をｘ（画素値ｘ）としたとき、画素値ｘの補正係数をｇ（ｘ）とすると、補正後の画素値ｆ（ｘ）は以下の［数１１］で表される。

　画素値の補正係数ｇ（ｘ）としては、例えば、画素値の測定値ｘの二次関数で近似することが考えられる。すなわち、基準となる画素値をｚとすると、画素値の補正係数ｇ（ｘ）は以下の［数１２］で表される。[数１２]より、基準となる画素値ｚにおける補正係数ｇ（ｚ）は１となる。なお、二次関数で近似したときの係数ａおよび係数ｂは、グラフ６ｂの線量と画素値の関係から求めることができる。

　なお、第１実施形態では画素値の補正係数ｇ（ｘ）を二次関数で近似する例を示したが、本発明はこのような形態に限定されない。例えば、画素値の補正係数ｇ（ｘ）を指数関数または対数などを用いて近似してもよいし、補正後の画素値ｆ（ｘ）を三次関数で近似するなどしてもよい。すなわち、補正後の画素値ｆ（ｘ）や画素値の補正係数ｇ（ｘ）を任意の関数で近似してもよい。また、線量と測定値の関係を複数点測定し、それらを補間して求めてもよい。

　また、補正後の画素値ｆ（ｘ）や補正係数ｇ（ｘ）を近似する関数に関して、二次元検出器１０６の画素毎に関数を変えてもよいし、二次元検出器１０６を複数の領域に分割して領域ごとで共通の関数を用いるようにしてもよい。ただし、計算を簡略化する観点から、二次元検出器１０６の全ての画素で共通の関数を用いるのが好ましい。また、画素値の測定値ｘを入力すると補正後の画素値ｆ（ｘ）を返すようなテーブルを用いて補正を行うように構成してもよい。

　なお、グラフ６ｂに示した線量と画素値の関係を測定する方法としては、Ｘ線発生装置１０１の管電流又は曝射時間を変更しながら撮影を行い、画素値を得る構成が用いられ得る。線量は管電流と曝射時間の積に比例するので、グラフ６ｂにおける理想の直線は、ある管電流と曝射時間で撮影したときの測定値と原点とを結ぶ直線となる。また、線量を変更する方法として、Ｘ線発生装置１０１から管電流一定で連続Ｘ線を曝射しながら、二次元検出器１０６の撮影間隔（露出時間）を変更して撮影する構成も好適に用いられる。この場合、線量は管電流と撮影間隔の積に比例するものとして、理想の直線を算出する。

　補正係数ｇ（ｘ）を求める方法は、グラフ６ｂに示した線量と画素値の関係を測定する方法に限定されない。例えば、グラフ６ａに示すような、被写体を置いたときの減弱率と線量の関係を測定し、減弱率が線量によらず一定になるような補正係数ｇ（ｘ）を定義するようにしてもよい。

　次に、エネルギーサブトラクション処理における線量依存補正の実行タイミングを説明する。以下では、特に、図３や図４で示したように、サンプルホールドを用いて撮影を行う構成へ上述した本実施形態の線量依存補正を適用する場合を説明する。すなわち、１ショットの放射線の曝射の間に複数回のサンプルホールドを行って得られた複数枚の放射線画像を取得して、エネルギーサブトラクション処理に提供する複数の画像を生成する構成における、線量依存補正の実行のタイミングについて説明する。

　図７Ａに、第１実施形態に係る補正処理の動作例を示す。まず、図５Ａに示したエネルギーサブトラクション処理における補正処理と同様に、オフセット補正と色補正を行い、被写体がない場合の低エネルギー画像Ｗ＿Ｌｏｗと被写体がない場合の高エネルギー画像Ｗ＿Ｈｉｇｈを得る。同様に、被写体がある場合の低エネルギー画像Ｘ＿Ｌｏｗと被写体がない場合の高エネルギー画像Ｘ＿Ｈｉｇｈを得る。撮像制御装置１０３は、オフセット補正と色補正を経て得られたこれらの画像に対して、グラフ６ｃで示した線量依存補正を行う。こうして、線量依存補正後の、被写体がない場合の低エネルギー画像Ｗ＿Ｌｏｗ'と被写体がない場合の高エネルギー画像Ｗ＿Ｈｉｇｈ'、被写体がある場合の低エネルギー画像Ｘ＿Ｌｏｗ'と被写体がある場合の高エネルギー画像Ｘ＿Ｈｉｇｈ'が得られる。

　その後、線量依存補正後の画像を用いて図５Ａで説明したゲイン補正を行うことで、低エネルギーにおける減弱率の画像Ｌと、高エネルギーにおける減弱率の画像Ｈが得られる。色補正後の画像に対して線量依存補正を行うことで、低エネルギーにおける減弱率Ｌと、高エネルギーにおける減弱率Ｈが線量に依存しなくなる（または依存性が低減する）。そのため、図５Ｂに示した信号処理によって算出される骨の厚みＢや軟部組織の厚みＳの推定値が、真値に近づくことが期待される。

　なお、図７Ａでは、色補正の実行後に線量依存補正を行っている（オフセット補正→色補正→線量依存補正→ゲイン補正の順番で処理が行われている）が、これに限られるものではない。例えば、図７Ｂに示されるように、色補正の実行前に線量依存補正を行う（オフセット補正→線量依存補正→色補正→ゲイン補正の順番で処理を行う）ようにしてもよい。

　図７Ｂにおいて、ＷＦ＿ＯＤＤは高エネルギーのＸ線４０２に対応する画像であり、ＷＦ＿ＥＶＥＮは低エネルギーのＸ線４０１、高エネルギーのＸ線４０２、低エネルギーのＸ線４０３の和に対応する画像である。したがって、低エネルギーのＸ線に対応する画像４０５を得るには、ＷＦ＿ＥＶＥＮからＷＦ＿ＯＤＤを減算すること、すなわち色補正を行う必要がある。図７Ｂでは、色補正に先立って線量依存補正が行われる。ＷＦ＿ＥＶＥＮとＷＦ＿ＯＤＤには線量依存があるため、図７Ｂに示されるように線量依存補正を行ってから色補正を行うことも好ましい処理順である。このとき、ＷＦ＿ＥＶＥＮとＷＦ＿ＯＤＤ、ＸＦ＿ＥＶＥＮとＸＦ＿ＯＤＤに対する線量依存補正の関数と係数を変えてもよい。ただし、計算を簡略化するという観点からは、ＷＦ＿ＥＶＥＮ、ＷＦ＿ＯＤＤ、ＸＦ＿ＥＶＥＮ、ＸＦ＿ＯＤＤで共通の関数と係数を用いて線量依存補正を行う構成が好ましい。

　＜第２実施形態＞
　第１実施形態では、減弱率の線量依存性を補正する線量依存補正を説明した。第２実施形態では、減弱率の厚みに依存する誤差およびスペクトルに依存する誤差に対する補正について補正説明する。なお、Ｘ線撮像システム及びそれを構成する各装置の構成は第１実施形態と同様である。

　図８のグラフ８ａに、減弱率の厚み依存の例を示す。第１実施形態（図６）に示した線量依存補正を行うことで、線量が変化しても減弱率が一定に保たれるようになる。しかしながら、管電圧やフィルタ等を変えず、すなわちスペクトルを変えずに被写体の厚みを変化させながら減弱率を測定したところ、グラフ８ａのように、減弱率の測定値が理想の曲線または直線から外れることが判明した。このような誤差の原因として、［数３］で示した減弱率の計算を行う際に用いるスペクトルＮ（Ｅ）や、エネルギーＥにおける骨の線減弱係数μ_Ｂ（Ｅ）、エネルギーＥにおける軟部組織の線減弱係数μ_Ｓ（Ｅ）が実際の値からずれていることが挙げられる。

　このような誤差を防ぐために、スペクトルＮ（Ｅ）はスペクトロメータなどを用いて実測した値を用いることが望ましい。しかしながら、Ｘ線撮像装置１０４の外装にＸ線が吸収されてスペクトルＮ（Ｅ）が変化したり、蛍光体１０５で一部のＸ線が透過することでスペクトルＮ（Ｅ）が変化したりするなどの課題がある。したがって、誤差を防ぐために、Ｘ線撮像装置１０４におけるＸ線の吸収や透過を考慮したスペクトルＮ（Ｅ）を用いることが望ましい。本願発明者が検討を行ったところ、以上に示したスペクトルＮ（Ｅ）の修正を行っても依然として許容できない誤差が存在する場合があることが判明した。

　そこで、第２実施形態の補正部１３２は、被写体の厚みに対する減弱率の変化が上記の所定の曲線または直線に一致するように減弱率の画像における各画素の減弱率を補正することにより、厚み依存性の補正を行う。例えば、図８のグラフ８ｂに、第２実施形態に係る厚み依存の補正係数を示す。第２実施形態では、減弱率の測定値に対してグラフ８ｂのような補正係数をかけることで、グラフ８ａにおける減弱率の測定値が理想の曲線または直線に近づくように補正される。例えば、減弱率の測定値をｙとしたとき、減弱率の補正係数をｉ（ｙ）とすると、補正後の減弱率ｈ（ｙ）は以下の［数１３］で表される。

　減弱率の補正係数ｉ（ｙ）としては、減弱率の測定値ｙの二次関数で近似する構成が好適に用いられるが、これに限られるものではない。例えば、線量依存補正と同様、補正後の減弱率ｈ（ｙ）や減弱率の補正係数ｉ（ｙ）を任意の関数で近似してもよい。また、二次元検出器１０６の画素毎に関数を変えてもよいし、二次元検出器１０６を複数の領域に分けて領域ごとで共通の関数を用いるようにしてもよい。ただし、計算を簡略化するためには、二次元検出器１０６の全ての画素で共通の関数を用いる構成が好ましい。

　次に、減弱率のスペクトル依存性に対する補正について説明する。図９のグラフ９ａに、減弱率のスペクトル依存を示す。第１実施形態（図６）で説明した線量依存補正を行うことで、線量が変化しても減弱率が一定に保たれるようになる。しかしながら、被写体の厚みを変えず、Ｘ線の管電圧を変えることでスペクトルＮ（Ｅ）の平均エネルギーを変化させながら減弱率を測定したところ、グラフ９ａのように、減弱率の測定値が理想の曲線または直線から外れることが判明した。このような誤差の原因は、グラフ８ａに関して上述した通りである。

　そこで、第２実施形態の補正部１３２は、放射線の平均エネルギーに対する減弱率の変化が所定の曲線または直線に一致するように、減弱率の画像における各画素の減弱率を補正するエネルギー依存性の補正を行う。図９のグラフ９ｂに、第２実施形態に係るスペクトル依存の補正係数を示す。第２実施形態では、減弱率の測定値に対してグラフ９ｂのような補正係数をかけることで、グラフ９ａにおける減弱率の測定値が理想の曲線または直線に近づくように補正を行う。例えば、減弱率の測定値をｙ、スペクトルＮ（Ｅ）の平均エネルギーをｅとしたとき、減弱率の補正係数をｋ（ｅ）とすると、補正後の減弱率ｊ（ｅ）は以下の［数１４］で表される。

　減弱率の補正係数ｋ（ｅ）は線量依存補正と同様、任意の関数で近似してよい。また、二次元検出器１０６の画素毎に関数を変えてもよいし、二次元検出器１０６の領域ごとで共通の関数を用いるようにしてもよい。ただし、計算を簡略化するためには、二次元検出器１０６の全ての画素で共通の関数を用いる構成が好ましい。

　以上、第２実施形態では、説明を簡略化するために、厚み依存性とスペクトル依存性を独立に補正する方法を示した。しかしながらこのような形態に限られるものではない。例えば、骨の厚みＢと軟部組織の厚みＳとスペクトルＮ（Ｅ）ごとに減弱率を実測し、理想の減弱率の曲線に一致するような補正係数ｌ（Ｂ，Ｓ，Ｎ（Ｅ））を求め、補正係数ｌ（Ｂ，Ｓ，Ｎ（Ｅ））を減弱率の測定値に適用して補正するようにしてもよい。ここで、理想の減弱率の曲線は、例えば、物質特性画像において分離または識別される被写体（上記では骨、軟部組織）について、当該被写体の減弱係数と厚みとＸ線エネルギーに対応したＸ線スペクトルとから算出される減弱率の曲線である。減弱係数、Ｘ線スペクトルは公知のデータベースから取得可能である。なお、このとき、全ての骨の厚みＢと軟部組織の厚みＳとスペクトルＮ（Ｅ）ごとに補正係数ｌ（Ｂ，Ｓ，Ｎ（Ｅ））を求めると、データ量が膨大となる。従って、複数の条件で取得した補正係数から補間によってその間の条件における補正係数を求めるようにしてもよい。以上に示した処理を、減弱率の補正と呼ぶものとする。

　また、上記では、減弱率の測定値に対してグラフ８ｂのような補正係数をかけることで、グラフ８ａにおける減弱率の測定値が理想の曲線または直線に近づくように補正を行う方法を示した。しかしながら本発明はこのような形態に限定されない。例えば、減弱率の測定値と計算値が一致するようにスペクトルＮ（Ｅ）を変形するようにしてもよい。すなわち、あるスペクトルＮ（Ｅ）において、複数種類の材質について、材質ごとに厚みの異なる複数の被写体を撮影し、測定された減弱率と［数２］または［数３］によって算出される減弱率が一致するように、放射線のスペクトルを変形するようにしてもよい。或いは、所定の材質、所定の厚みの被写体を異なる複数のＸ線エネルギーで撮影し、測定された減弱率と［数２］または［数３］によって算出される減弱率が一致するように、放射線のスペクトルを変形するようにしてもよい。これらの場合、変形された放射線スペクトルは撮像制御装置１０３により保持され、信号処理で用いられることになる。したがって、上述したようなゲイン補正後の画像に対して補正を行う処理は不要となり、第１実施形態で説明した補正処理（図７Ａ、図７Ｂ）を行った後に図５Ｂに示される信号処理を実行すればよく、その信号処理において、変形されたスペクトルＮ（Ｅ）が用いられることになる。なお、放射線スペクトルを変形するための減弱率の測定において、第１実施形態で説明した、線量依存の誤差を補正する補正処理（図７Ａ、図７Ｂ）を行って減弱率画像を得ることが好ましい。

　図１０に、第２実施形態に係る撮像制御装置１０３が行う補正処理の動作例を示す。図７Ｂで示したように、オフセット補正→線量依存補正→色補正→ゲイン補正の順番で補正が行われる。なお、図７Ａに示した順番（オフセット補正→色補正→線量依存補正→ゲイン補正）が用いられてもよい。図１０ではさらに、ゲイン補正の後に減弱率の補正（減弱率の厚み依存および減弱率のスペクトル依存の補正）が行われ減弱率の画像Ｈ，Ｌがそれぞれ減弱率の画像Ｈ'、Ｌ'となる。減弱率の補正を行うことで、低エネルギーにおける減弱率Ｌと、高エネルギーにおける減弱率Ｈが厚みやスペクトルに依存しなくなる（または、依存性が低減する）ため、図５Ｂに示した信号処理によって算出した骨の厚みＢや軟部組織の厚みＳの推定値が、真値に近づくことが期待される。

　＜第３実施形態＞
　第１実施形態、第２実施形態では、センサ出力やスペクトルを補正する方法について説明した。第３実施形態、第４実施形態では、Ｘ線発生装置１０１とＸ線撮像装置１０４の間の被写体以外の構成物による放射線（Ｘ線）の吸収が、構成物の個体ばらつきにより生じる誤差を低減する校正方法を提案する。Ｘ線吸収に影響する構成物の例として、第３実施形態では、Ｘ線を可視光に変換する蛍光体１０５を例示して説明する。第３実施形態に係る撮像システムの構成および画素２０の等価回路図は、第１実施形態（図１、図２）と同様である。

　図１１の構成１１ａは、第３実施形態におけるＸ線撮影における構成を示す。Ｘ線発生装置１０１とＸ線撮像装置１０４との間に被写体１１０１が配置された状態でＸ線発生装置１０１からＸ線を照射することで、被写体１１０のＸ線画像が得られる。Ｘ線撮像装置１０４には、Ｘ線を可視光に変換する蛍光体１０５が配されており、蛍光体１０５によるＸ線の吸収が、エネルギーサブトラクション処理の精度に影響を及ぼす。Ｘ線発生装置１０１から曝射されたＸ線のエネルギーＥにおけるスペクトルをＮ_０（Ｅ）、蛍光体１０５のエネルギーＥにおける蛍光体の線減弱係数をμ_Ｃ（Ｅ）、厚さをｄ_Ｃ、充填率をＰ_Ｃとすると、以下の［数１５］に示される式が成り立つ。

　［数１５］のＸ線スペクトルＮ（Ｅ）は、蛍光体１０５のＸ線吸収を考慮したスペクトルであり、Ｘ線スペクトルＮ_０（Ｅ）と蛍光体１０５のＸ線吸収に関わるパラメータから得られる。Ｘ線発生装置１０１から照射されるＸ線のＸ線スペクトルＮ_０（Ｅ）は、シミュレーション又は実測により得られる。また、蛍光体１０５のエネルギーＥにおける線減弱係数μ_Ｃ（Ｅ）、厚さｄ_Ｃ、充填率Ｐ_Ｃは設計値から得られる。さらに、エネルギーＥにおける任意の物質の線減弱係数μ（Ｅ）は、ＮＩＳＴなどのデータベースから得られる。したがって、任意の物質の厚さＤ、Ｘ線のスペクトルＮ（Ｅ）における減弱率Ｉ／Ｉ_０を計算することが可能である。なお、［数１５］については、簡単のため線源弱係数のみを用いているが、実際の計算ではエネルギー吸収係数を考慮した方がより精度の高い結果が得られる。蛍光体１０５を透過しなかったＸ線（［数１５］によって示されるＸ線）がすべてセンサ出力に変換されるわけではない。したがって、センサ出力に変換されないＸ線を考慮すること、すなわちエネルギー吸収係数を考慮することで、精度向上が期待される。

　図１１の構成ｂは、第３実施形態における補正係数を取得するための構成を示す。Ｘ線発生装置１０１とＸ線撮像装置１０４の間に校正用の試料１１０２が配置されている。このとき、Ｘ線の線減弱係数、密度、厚さが既知の試料１１０２を用いれば、［数３］をもとに試料１１０２の減弱率Ｉ／Ｉ_０の計算値が得られる。また、構成１１ｂで実際にＸ線撮影を行い、Ｘ線画像を取得することで試料１１０２の減弱率Ｉ／Ｉ_０の実測値が得られる。

　図１１の１１ｃに減弱率Ｉ／Ｉ_０の実測の例を示す。取得されたＸ線画像１１２１は画素２０の数に応じた、Ｘ列×Ｙ行の２次元情報である。Ｘ線画像１１２１内の領域１１２２はＸ線が試料１１０２を透過してＸ線撮像装置１０４に入射した領域１１２４の一部であり、領域１１２３はＸ線が試料１１０２を透過せずにＸ線撮像装置１０４に入射した領域１１２５の一部である。このとき、領域１１２２内の画素２０の出力平均値Ｉ、領域１１２３内の画素２０の出力平均値Ｉ_０を求めることで、減弱率Ｉ／Ｉ_０の実測値が得られる。なお、減弱率Ｉ／Ｉ_０を得るまでに、第１実施形態で説明したようなオフセット補正、色補正、線量依存補正、第２実施形態で説明したような厚み依存補正、エネルギー依存補正が行われてもよい。

　以上のようにして試料１１０２について得られた減弱率Ｉ／Ｉ_０の計算値と実測値は、一般には一致しない。減弱率の計算値に、蛍光体１０５の個体ばらつきによる誤差要因が含まれるためである。そこで、試料１１０２の減弱率Ｉ／Ｉ_０の計算値と実測値が一致するような、蛍光体１０５の線減弱係数μ_Ｃ'（Ｅ）、厚さｄ_Ｃ'、充填率Ｐ_Ｃ'を求める。線減弱係数μ_Ｃ'（Ｅ）、厚さｄ_Ｃ'、充填率Ｐ_Ｃ'は、例えば、試料１１０２の減弱率Ｉ／Ｉ_０の計算値と実測値の差分をもとに求めてもよいし、反復法で計算値を実測値に収束させて求めてもよい。

　図１２に、第３実施形態による補正処理のフローチャートを示す。Ｓ１２０１において、取得部１３１は、試料１１０２のＸ線画像として、図１１の１１ｃに示したようなＸ線画像１１２１を取得する。Ｓ１２０２において、補正部１３２は、Ｘ線画像１１２１の領域１１２２と領域１１２３の画素値を用いて減弱率（Ｉ／Ｉ_０）を求め、実測値として保持する。一方、補正部１３２は、Ｓ１２０３において［数１５］で使用されるパラメータ（線減弱係数μ_Ｃ（Ｅ）、厚さｄ_Ｃ、充填率Ｐ_Ｃ）の設計値を取得する。そして、補正部１３２は、Ｓ１２０４において、Ｓ１２０３で取得したパラメータと［数１５］を用いて減弱率（Ｉ／Ｉ_０）を計算する。

　Ｓ１２０５において、補正部１３２は、Ｓ１２０２で取得された減弱率の実測値とＳ１２０４で取得された減弱率の計算値との差を計算する。そして、Ｓ１２０６において、補正部１３２は、Ｓ１２０５で計算された差に基づいて、対象物質（本例では蛍光体１０５）のパラメータの設計値に対する誤差を計算する。Ｓ１２０７において、補正部１３２は、Ｓ１２０６で算出された誤差を蛍光体１０５のパラメータに反映することにより、補正された線減弱係数μ_Ｃ'（Ｅ）、厚さｄ_Ｃ'、充填率Ｐ_Ｃ'を得る。蛍光体１０５のこれらの補正されたパラメータは、設計値に対する個体ばらつきが補正されたものである。従って、これら補正されたパラメータを［数１５］に代入することで、蛍光体１０５の固体ばらつきが補正されたスペクトルＮ（Ｅ）が得られる。こうして得られたスペクトルＮ（Ｅ）を用いて［数４］を解くことで、算出した骨の厚みＢや軟部組織の厚みＳの推定値が、真値に近づくことが期待される。

　Ｓ１２０８において、取得部１３１は、複数の異なるＸ線エネルギーにより被写体を撮影した複数のＸ線画像を取得し、補正部１３２は、例えば図７Ａ、（ｂ）で説明した処理を行って複数の異なるＸ線エネルギーに対応した複数の減弱率画像を得る。Ｓ１２０９において、信号処理部１３３は、Ｓ１２０７で生成された補正されたパラメータと［数１５］を用いて、補正後のＸ線スペクトルＮ（Ｅ）を生成し、補正後のＸ線スペクトルＮ（Ｅ）を用いてエネルギーサブトラクション処理を行う。こうして、補正されたＸ線スペクトルを用いたエネルギーサブトラクション処理により物質分離演算が行われ、例えば、骨と軟部組織の厚み画像が得られる。

　なお、減弱率Ｉ／Ｉ_０の実測値（および計算値）のデータは、必要に応じてパターンを変えて複数取得するようにしてもよい。パターンを変える例としては、Ｘ線発生装置１０１で照射するＸ線のエネルギーを変える、試料１１０２の材料を変える、試料１１０２の厚みを変えること等があげられる。また、第３実施形態では、蛍光体１０５のＸ線吸収にかかわる３つのパラメータ（μ_Ｃ（Ｅ）、ｄ_Ｃ、Ｐ_Ｃ）を補正する場合について説明したが、これに限定されるものではない。例えば、補正の対象となるパラメータはそれら３つのパラメータのうちの１つまたは２つであってもよいし、補正の対象となるパラメータをユーザが任意に選択できるようにしてもよい。例えば、蛍光体の厚み（ｄ_Ｃ）のみを補正するようにしてもよい。補正対象となるパラメータが増える程、必要な減弱率Ｉ／Ｉ_０の実測値（および計算値）の数が増加するため、物質の設計ばらつきが小さく減弱率への影響が小さいパラメータは補正対象から除くようにしてもよい。また、質量減弱係数、密度といったＸ線の減弱に影響する他のパラメータを補正してもよい。

　なお、設計値（μ_Ｃ（Ｅ）、ｄ_Ｃ、Ｐ_Ｃ）と補正後の値（μ_Ｃ'（Ｅ）、ｄ_Ｃ'、Ｐ_Ｃ'）の差分Δμ_Ｃ（Ｅ）、厚さΔｄ_Ｃ、充填率ΔＰ_Ｃを補正値としたとき、それら補正値の絶対値（補正範囲）は設計ばらつきの範囲に制限するのが好ましい。補正値が非現実的な値をとることが防止されるからである。

　また、第３実施形態では、試料１１０２の減弱率Ｉ／Ｉ_０の実測値を取得する際に、Ｘ線画像１１２１に試料ありの領域１１２４、試料なしの領域１１２５を設け、各領域内において平均値を算出する方法を用いたが、これに限定されるものではない。例えば、試料ありの画像、試料なしの画像をそれぞれ取得するように撮影を複数回に分けて行うようにしてもよい。また、このとき、補正値を画素２０ごとに取得し、画素単位で上述したパラメータの補正を行うようにしてもよい。このとき、補正値はＸ列×Ｙ行の２次元情報となり、対象物体（本例では蛍光体１０５）の面内のばらつきも補正できる。

　＜第４実施形態＞
　第３実施形態では、蛍光体によるＸ線の吸収を考慮する構成を説明した。第４実施形態では、さらに、撮影の際に、Ｘ線発生装置１０１とＸ線撮像装置１０４（蛍光体１０５）の間に、被写体以外にＸ線を吸収する構成として、付加フィルタおよび／または散乱線除去用のグリッドが存在する場合について説明する。第４実施形態に係る撮像システムのブロック図および、画素２０の等価回路図は、第１実施形態（図１、図２）と同様である。

　図１３の１３ａに、第４実施形態におけるＸ線撮影の構成例(構成１３ａ)を示す。Ｘ線発生装置１０１とＸ線撮像装置１０４の間に被写体１１０１が配置される。また、Ｘ線発生装置１０１には付加フィルタ１３０１が装着されており、Ｘ線撮像装置１０４には散乱線除去用のグリッド１３０２が装着されている。付加フィルタ１３０１には、ＡｌやＣｕといった金属性フィルタが好適に用いられ、低エネルギーのＸ線領域をカットする。グリッド１３０２は被写体１１０１で発生し得る散乱線を除去する。Ｘ線発生装置１０１から曝射されたＸ線のスペクトルをＮ_０（Ｅ）とし、付加フィルタ１３０１とグリッド１３０２を透過した後のＸ線のスペクトルをＮ'（Ｅ）とする。また、フィルタの厚さをｄ_Ｆ、グリッドの厚さをｄ_Ｇ、エネルギーＥにおける骨の線減弱係数をμ_Ｆ（Ｅ）、エネルギーＥにおけるグリッドの線減弱係数をμ_Ｇ（Ｅ）とすると、蛍光体１０５に到達するＸ線に関して、以下の［数１６］の式が成り立つ。

　さらに、蛍光体１０５の厚さをｄ_Ｃ、エネルギーＥにおける蛍光体の線減弱係数をμ_Ｃ（Ｅ）、充填率をＰ_Ｃとすると、以下の［数１７］の式が成り立つ。［数１７］のＸ線スペクトルＮ（Ｅ）は、付加フィルタ１３０１、グリッド１３０２、蛍光体１０５のＸ線吸収を考慮したスペクトルである。

　[数１６]、［数１７］により得られたＸ線スペクトルＮ（Ｅ）を用いて、第３実施形態と同様に［数４］の非線形連立方程式を解くことで骨の厚みＢや軟部組織の厚みＳが得られる。図１３の構成１３ａにおいて［数４］を解くことにより得られる骨や軟部組織の厚みの推定値は、真値との間に誤差が生じる。第３実施形態で説明したような蛍光体１０５の個体ばらつきに加え、付加フィルタ１３０１、グリッド１３０２の個体ばらつきが含まれるためである。

　図１３の１３ｂに、第４実施形態における補正係数取得の構成を示す。Ｘ線発生装置１０１とＸ線撮像装置１０４の間に試料１１０２が配置される。試料１１０２には、Ｘ線の源弱係数や密度、厚さが既知の物質が用いられる。補正処理の流れは第３実施形態（図１２）と同様である。補正部１３２は、付加フィルタ１３０１、グリッド１３０２、蛍光体１０５および試料１１０２の設計値（パラメータ）を用いて［数４］、［数１６］、[数１７]をもとに試料１１０２の減弱率Ｉ／Ｉ_０の計算値を得る（Ｓ１２０３、Ｓ１２０４）。また、補正部１３２は、図１３の構成１３ｂによりＸ線撮影を行い、Ｘ線画像を取得することで試料１１０２の減弱率Ｉ／Ｉ_０の実測値を得る（Ｓ１２０１、Ｓ１２０２）。

　補正部１３２は、第３実施形態と同様の方法で、補正後のパラメータを求める（Ｓ１２０５～Ｓ１２０７）。補正後のパラメータは、例えば、蛍光体１０５の線減弱係数μ_Ｃ'（Ｅ）、厚さｄ_Ｃ'、充填率Ｐ_Ｃ'、付加フィルタ１３０１の線減弱係数μ_Ｆ'（Ｅ）、厚さｄ_Ｆ'、グリッド１３０２の線減弱係数μ_Ｇ'（Ｅ）、厚さｄ_Ｇ'である。信号処理部１３３は、これら補正後のパラメータを［数１６］、［数１７］に代入して得られたスペクトルＮ（Ｅ）を用いて［数４］を解くことで、骨の厚みＢ、軟部組織の厚みＳの推定値を得る（Ｓ１２０８、Ｓ１２０９）。このように、補正されたＸ線スペクトルＮ（Ｅ）を用いてエネルギーサブトラクション処理を行うことにより、算出された骨の厚みＢや軟部組織の厚みＳの推定値が、真値に近づくことが期待される。

　以上、第４実施形態ではＸ線撮影の付属品として付加フィルタとグリッドを装着した場合について説明したが、本発明はこのような形態に限定されない。例えば、付加フィルタとグリッドのいずれか一方のみが装着されてもよいし、さらに別の付属部品が追加されてもよい。また、被写体を除き、Ｘ線発生装置１０１と蛍光体１０５の間にＸ線を吸収し得る物体がある場合は、その物体を補正対象としてもよい。そのような物体としては、例えば、Ｘ線撮像装置の外装部品などがあげられる。

　＜第５実施形態＞
　第３、第４実施形態では、センサ出力である減弱率Ｉ／Ｉ_０をもとに、被写体以外の構成（物体）のＸ線吸収に関わるパラメータを校正し、Ｘ線スペクトルを補正した。第５実施形態では、複数のＸ線エネルギーに対応した複数の減弱率画像を用いたエネルギーサブトラクション処理の処理結果をもとに、上記校正を行う。第５実施形態に係る撮像システムのブロック図および、画素２０の等価回路図は、第１実施形態（図１、図２）と同様である。

　第５実施形態による補正係数取得のための構成は、第３実施形態（構成１１ｂ）と同様である。すなわち、Ｘ線発生装置１０１とＸ線撮像装置１０４の間に校正用の試料１１０２が配置される。ここで、試料１１０２は、骨と軟部組織の部分（または、それらと同等のＸ線減弱係数を有する材料）のみで構成されており、各部分の厚みは既知である。図１１の構成１１ｂでＸ線撮影を行い、Ｘ線画像を取得することで校正用の試料１１０２の減弱率Ｉ／Ｉ_０の実測値が得られる。この状態で、高エネルギーにおける減弱率画像（高エネルギー画像）と低エネルギーにおける減弱率画像（低エネルギー画像）を得ることで、試料１１０２における骨の厚さの画像Ｂと軟部組織の厚さの画像Ｓを得ることができる。こうして得られた画像Ｂ、画像Ｓと試料１１０２の実際の各部分の厚みは、Ｘ線スペクトルＮ_Ｈ（Ｅ）、Ｎ_Ｌ（Ｅ）の誤差によって一致しない。そのため、Ｘ線スペクトルの校正を行う必要がある。第３、第４実施形態の校正方法では、減弱率Ｉ／Ｉ_０の実測値とＸ線スペクトルから計算される値が一致するようにＸ線スペクトルを変形した。しかしながら、減弱率の誤差が小さくなっても厚み誤差（エネルギーサブトラクション処理の結果の誤差）が小さくならない場合がある。

　第５実施形態では、図１４に示すような処理でＸ線スペクトルの校正を行う。図１４は、第５実施形態による補正処理を示すフローチャートである。Ｓ１４０１，Ｓ１４０２において、取得部１３１は、図１１の構成１１ｂにおいて複数の異なるＸ線エネルギーでＸ線撮影された複数のＸ線画像を取得する。補正部１３２は、取得部１３１が取得した複数の異なるＸ線エネルギーに対応した複数のＸ線画像から、高エネルギーの減弱率画像（高エネルギー画像）と低エネルギーの減弱率画像（低エネルギー画像）を生成する。Ｓ１４０３において、処理部１３３は、Ｓ１４０１とＳ１４０２で取得された高エネルギー画像と低エネルギー画像を用いてエネルギーサブトラクション処理を行い、試料１１０２の骨の厚みＢと軟部組織の厚みＳを計算する。Ｓ１４０４において、補正部１３２は、厚みＢおよび厚みＳの計算値と、試料１１０２の厚みＢおよび厚みＳの実測値との誤差を計算する。そして、Ｓ１４０５において、Ｓ１４０４で計算される誤差が小さくなるようにＸ線スペクトルＮ_Ｈ（Ｅ）、Ｎ_Ｌ（Ｅ）を補正する。

　次に、Ｓ１４０６において、取得部１３１、補正部１３２は、被写体１１０１のＸ線画像を撮影して異なる複数のエネルギーに対応した複数のＸ線画像を取得し、これらに基づいて複数の異なるエネルギーに対応した複数の減弱率画像を取得する。Ｓ１４０７において、処理部１３３は、補正されたスペクトルＮ_Ｈ（Ｅ）、Ｎ_Ｌ（Ｅ）を用いて物質分離を行う。以上の処理によって、被写体１１０１について、より真値に近い分離結果が得られると考えられる。

　Ｓ１４０５におけるＸ線スペクトルの補正方法の一例を示す。Ｘ線源であるＸ線発生装置１０１とセンサであるＸ線撮像装置１０４との間に、減弱係数μ（Ｅ）、厚みｄの想定していない物体が挟まっていると仮定する。この場合、元のＸ線スペクトルをＮ（Ｅ）とすると、校正後のＸ線スペクトルＮｃ（Ｅ）との間に以下の［数１８］が成り立つ。

　Ｎｃ（Ｅ）を用いて試料１１０２の厚みＢ、厚みＳを計算し、それぞれの計算値と実測値との誤差が小さくなるようなμ（Ｅ）、ｄを求める。このような処理を、
　第５実施形態では、Ｘ線源とセンサ間に減弱係数と厚みが未知の物質が挟まっているという物理モデルをもとにスペクトルの校正を行った。しかしながら本発明はこのような形態に限定されない。ランダムにスペクトルを変形しても良いし、実効原子番号及び面密度の定まった特定の物質が挟まっていると仮定しても良い。

　スペクトル校正における最適化の手順ついて説明する。まず、実効原子番号Ｚ、面密度Ｄの物質の値をランダムに決定する。この時、Ｚ、Ｄは現実的な値をとるように制限するのがよい。次に、決定したＺ、Ｄをもとに、高、低エネルギーのスペクトルを式（１７）に従って変形する。次に、変形したスペクトルをもとにＺ、Ｄを計算し、実測とのＲＭＳＥ（Root Mean Square Error）を求める。以上のようなスペクトルの変形を既定の回数繰り返し、最もＲＭＳが小さいスペクトルを求める。以上を既定の回数繰り返すことで、複数のスペクトルを求め、それらのスペクトルを平均する。平均されたスペクトルを校正前スペクトルとして、同様の手順を既定の回数繰り返す。以上の処理により、スペクトルの最適化を行う。なお、本発明は上述の形態に限定されるものではない。例えば、原子番号Ｚおよび厚みＤは必ずしも現実的な値をとる必要はないし、必ずしもスペクトルの変形が繰り返されなくてもよい。

　第５実施形態におけるスペクトルＮ（Ｅ）は［数１５］に示す通り、蛍光体１０５のＸ線吸収を考慮したスペクトルを示しており、撮影を行った際の構成をもとに、Ｘ線発生装置１０１から発生するＸ線のスペクトルを変形させた値を用いている。本実施形態ではスペクトルＮ（Ｅ）を変形する場合について述べたが、本発明はこのような形態に限定されず、スペクトル変形に使用したパラメータを補正しても良い。例えば、蛍光体の厚みｄ_Ｃや充填率Ｐ_Ｃを補正するようにしても良い。

　＜第６実施形態＞
　第５実施形態においては、一個のＸ線フォトンが蛍光体に吸収されたときに生じるセンサ出力は、Ｘ線フォトンのエネルギーＥに比例するものとしていた。すなわち、Ｘ線からセンサ出力への変換効率にはエネルギー依存がないと仮定されていた。第６実施形態では、センサの変換効率のエネルギー依存を考慮した場合の補正について説明する。第６実施形態に係る撮像システムのブロック図および、画素２０の等価回路図は、第１実施形態（図１、図２）と同様である。

　Ｘ線フォトンのエネルギーをＥ、Ｘ線のスペクトルをＮ（Ｅ）、センサ出力の変換効率をＣ（Ｅ）、任意物質の厚さをＤ、任意物質の線減弱係数をμ（Ｅ）、減弱率をＩ／Ｉ_０とすると、以下の［数１９］に示される式が成り立つ。

　［数１９］で、Ｎ_１（Ｅ）＝Ｃ（Ｅ）Ｎ（Ｅ）とすれば、［数３］と同じ形になる。従って、第５実施形態で用いられた試料１１０２をＸ線撮影して得られた高エネルギー画像と低エネルギー画像について第１実施形態と同様の計算（エネルギーサブトラクション処理）を行うことで、試料１１０２の厚みＢ、Ｓが得られる。この計算によって得られた試料１１０２の厚みＢ、Ｓと、試料１１０２の実際の厚みＢ、Ｓとの誤差をもとにして変換効率Ｃ（Ｅ）を校正できる。

　Ｃ（Ｅ）の校正方法について説明する。校正前の変換効率Ｃ（Ｅ）はＣ（Ｅ）＝１とみなすことができる。これを適当な関数でフィッティングすることで校正する。例えば、Ｃ（Ｅ）＝ａＥ＋ｂとする。変数ａ、ｂをランダムに決定して厚みＢ、Ｓを計算し、実際の厚みＢ、ＳとのＲＭＳＥ（Root Mean Square Error）を求める。以上を既定の回数繰り返してＲＭＳＥが最も小さいａ、ｂを求めることで校正された変換効率Ｃｃ（Ｅ）が得られる。処理部１３３は、被写体１１０１のＸ線画像を撮影して、校正された変換効率Ｃｃ（Ｅ）を用いて物質分離を行う。こうして、Ｃｃ（Ｅ）Ｎ（Ｅ）を校正されたＸ線スペクトルとして用いてエネルギーサブトラクション処理を行うことによって、被写体１１０１について、より真値に近い分離結果（厚みＢ、厚みＳ）が得られることが期待される。

　なお、第６実施形態ではＣ（Ｅ）を一次関数と仮定してＢ，ＳのＲＭＳＥをもとに校正を行ったが、本発明はこのような形態に限定されない。例えば、Ｃ（Ｅ）の関数は二次関数や自然対数などの線形和でもよい。また、Ｂ、ＳのＲＭＳＥではなく、実測値と径産地との間の差分や相対誤差をもとにＣ（Ｅ）の校正が行われても良い。また、第６実施形態では、Ｂ、Ｓの計算値と実測値の誤差をもとにＣ（Ｅ）を補正する場合について述べたが、本発明はこのような形態に限定されない。例えば、Ｂのみの厚みをもとに補正しても良いし、Ｃ（Ｅ）だけでなく、Ｎ（Ｅ）を同時に補正しても良い。

　なお、第１～第６実施形態では、エネルギーサブトラクション処理により骨の厚さＢと軟部組織の厚さＳを算出していたが、本発明はこのような形態に限定されない。例えば、水の厚さＷと造影剤の厚さＩが算出されてもよい。すなわち、本発明は、任意の二種類の物質の厚さへの分解に適用され得る。また、低エネルギーにおける減弱率の画像Ｌと高エネルギーにおける減弱率の画像Ｈから、エネルギーサブトラクション処理により実効原子番号Ｚの画像と面密度Ｄの画像を計算する構成でもよい。実効原子番号Ｚとは混合物の等価的な原子番号のことであり、面密度Ｄとは被写体の密度［ｇ／ｃｍ^３］と被写体の厚み［ｃｍ］の積である。さらに、画像処理部１３４は、実効原子番号Ｚと面密度Ｄを用いて仮想単色Ｘ線画像を生成してもよい。また、画像処理部１３４は、複数のエネルギーＥ_Ｖで生成した複数の仮想単色Ｘ線画像を合成することで、合成Ｘ線画像を生成してもよい。合成Ｘ線画像とは、任意のスペクトルのＸ線を照射した場合に得られることが想定される画像のことである。

　また、上記の第１～第６実施形態では、Ｘ線撮像装置１０４として蛍光体を用いた間接型のＸ線センサを用いた。しかしながら本発明はこのような形態に限定されない。例えばＣｄＴｅ等の直接変換材料を用いた直接型のＸ線センサを用いてもよい。すなわち、Ｘ線センサは、間接型・直接型のどちらでもよい。また、第１、第２実施形態では、例えば、図４の動作においてＸ線発生装置１０１の管電圧を変化させていた。しかしながら本発明はこのような形態に限定されない。Ｘ線発生装置１０１のフィルタを時間的に切り替えるなどして、Ｘ線撮像装置１０４に曝射されるＸ線のエネルギーを変化させるようにしてもよい。すなわち、Ｘ線撮像装置１０４へ曝射されるＸ線のエネルギーを変更させる方法は、どのような方法は、何等限定されない。また、第２実施形態では、第１実施形態による線量依存補正の実施が前提になっているが、これに限られるものではない。第２実施形態で説明した減弱率の補正のみを行うようにしてもよい。また、厚み依存性に対する補正とエネルギー依存性に対する補正の一方のみを実行する構成であってもよい。

　また、第１～第６実施形態では、Ｘ線のエネルギーは２つであったが、これに限定されるものではない。Ｘ線のエネルギーが３つ以上の場合についても、例えば、図１０に示したように、オフセット補正→線量依存補正→色補正→ゲイン補正→減弱率の補正の順番で補正を行う構成を適用できる。すなわち、３つ以上のＸ線エネルギーから得られるＸ線画像についても上記実施系他の処理が適用可能である。

　また、図１０では、図７Ｂによる補正処理順を適用しているが、図７Ａに示した補正処理順（オフセット補正→色補正→線量依存補正→ゲイン補正→減弱率の補正）が適用されてもよい。また、第１～第６実施形態では、Ｘ線のエネルギーを変化させることで異なるエネルギーの画像を得ていたが、本発明はこのような形態に限定されない。例えば、複数の蛍光体１０５および二次元検出器１０６を重ねることで、Ｘ線の入射方向に対して前面の二次元検出器と背面の二次元検出器から異なるエネルギーの画像を得る、積層型の構成としてもよい。この場合、補正部１３２による色補正は不要となる。

　また、第１～第６実施形態では、Ｘ線撮影システムの撮像制御装置１０３を用いてエネルギーサブトラクション処理を行っていた。しながらこの本発明はこのような形態に限定されない。例えば、撮像制御装置１０３で取得した画像を別のコンピュータに転送して、エネルギーサブトラクション処理を行ってもよい。例えば、取得した画像を医療用のＰＡＣＳを介して別のパソコンに転送し、エネルギーサブトラクション処理を行ってから表示する構成としてもよい。すなわち、上記実施形態で説明した補正処理を行う装置は、撮影装置とセットでなくてよい（画像ビューアでもよい）。

　＜他の実施形態＞
　本発明は、上述の実施形態の１以上の機能を実現するプログラムを、ネットワーク又は記憶媒体を介してシステム又は装置に供給し、そのシステム又は装置のコンピュータにおける１つ以上のプロセッサーがプログラムを読出し実行する処理でも実現可能である。また、１以上の機能を実現する回路（例えば、ＡＳＩＣ）によっても実現可能である。

　本発明は上記実施の形態に制限されるものではなく、本発明の精神及び範囲から離脱することなく、様々な変更及び変形が可能である。従って、本発明の範囲を公にするために、以下の請求項を添付する。

　本願は、２０１９年６月１１日提出の日本国特許出願特願２０１９－１０９０１９および２０２０年５月２８日提出の日本国特許出願特願２０２０－０９３６７３を基礎として優先権を主張するものであり、それらの記載内容の全てを、ここに援用する。

Claims (26)

1. 　被写体に放射線を照射して撮影を行うことで得られた互いに異なる複数の放射線エネルギーに対応した複数の放射線画像を取得する取得手段と、
   　前記取得手段により取得された前記複数の放射線画像からエネルギーサブトラクション処理で用いられる複数の減弱率画像を生成する生成手段と、
   　前記放射線の線量、前記被写体の厚み、および前記放射線のエネルギーの少なくともいずれかに依存して生じる減弱率の誤差を低減するように、放射線画像または減弱率画像を補正する補正手段と、
   　前記補正手段による補正後の前記複数の減弱率画像を用いて前記エネルギーサブトラクション処理により物質特性画像を生成する処理手段と、を備えることを特徴とする画像処理装置。
2. 　前記補正手段は、画素値と線量が比例するように補正係数を用いて放射線画像の画素値を補正することを特徴とする請求項１に記載の画像処理装置。
3. 　前記補正手段は、オフセット補正を実行した後の放射線画像の画素値を、前記補正係数を用いて補正することを特徴とする請求項２に記載の画像処理装置。
4. 　前記取得手段は、前記複数の放射線画像として、１ショットの放射線の曝射の間に複数回のサンプルホールドを行って得られた放射線画像を取得し、
   　前記補正手段は、前記複数の放射線画像の減算により、サンプルホールドによって区切られる特定の期間の放射線に対応した画像を得る色補正を含み、前記色補正の実行後の画像に前記補正係数を用いた補正を行うことを特徴とする請求項２または３に記載の画像処理装置。
5. 　画素値をｘ、前記補正係数をｇ（ｘ）とした場合に、前記補正手段はｆ（ｘ）＝ｘ×ｇ（ｘ）により補正後の画素値ｆ（ｘ）を取得し、補正係数ｇ（ｘ）が二次関数で表されることを特徴とする請求項２乃至４のいずれか１項に記載の画像処理装置。
6. 　前記補正手段は、画素値ｘと補正後の画素値ｆ（ｘ）の関係を格納しているテーブルを参照して前記補正を行うことを特徴とする請求項２乃至５のいずれか１項に記載の画像処理装置。
7. 　前記補正手段は、厚みに依存する誤差または放射線のエネルギーに依存する誤差を低減するように前記複数の減弱率画像の画素値を補正することを特徴とする請求項１乃至６のいずれか１項に記載の画像処理装置。
8. 　前記補正手段は、被写体の厚みに対する減弱率の変化が所定の曲線または直線に一致するように前記複数の減弱率画像の画素値を補正することを特徴とする請求項７に記載の画像処理装置。
9. 　前記補正手段は、減弱率の関数である補正係数を前記減弱率に乗じることで前記厚み依存性の補正を行うことを特徴とする請求項８に記載の画像処理装置。
10. 　前記補正手段は、放射線の平均エネルギーに対する減弱率の変化が所定の曲線または直線に一致するように、前記複数の減弱率画像の画素値を補正することを特徴とする請求項７乃至９のいずれか１項に記載の画像処理装置。
11. 　前記補正手段は、放射線の平均エネルギーの関数である補正係数を前記減弱率に乗じることで、前記減弱率のエネルギー依存性を補正することを特徴とする請求項１０に記載の画像処理装置。
12. 　エネルギーサブトラクション処理で用いられる放射線スペクトルを、計算により得られた画素値と実測により得られた画素値との差に基づいて補正する補正手段と、
    　被写体に放射線を照射して撮影を行うことで得られた互いに異なる複数の放射線エネルギーに対応した複数の放射線画像を取得する取得手段と、
    　前記取得手段により取得された複数の放射線画像から前記エネルギーサブトラクション処理に用いられる複数の減弱率画像を生成する生成手段と、
    　前記複数の減弱率画像と前記補正手段により補正された放射線スペクトルを用いて前記エネルギーサブトラクション処理により物質特性画像を生成する処理手段と、を備えることを特徴とする、画像処理装置。
13. 　前記補正手段は、複数の厚みの試料について取得される放射線画像から測定される減弱率および複数の放射線エネルギーで取得される試料の放射線画像から測定される減弱率の少なくとも何れかと、放射線スペクトルを用いて計算された前記試料の減弱率とが一致するように放射線スペクトルを補正することを特徴とする請求項１２に記載の画像処理装置。
14. 　前記補正手段は、前記放射線スペクトルを用いて計算された減弱率と実測により得られた減弱率、または、前記エネルギーサブトラクション処理により計算された物質特性と実測により得られた物質特性に基づいて、前記放射線スペクトルを補正する、ことを特徴とする請求項１２に記載の画像処理装置。
15. 　前記補正手段は、複数の放射線エネルギーのそれぞれの放射線スペクトルを、前記被写体以外の構成により生じる放射線の吸収に基づいて補正することを特徴とする請求項１４に記載の画像処理装置。
16. 　前記補正手段は、線減弱係数および厚さが既知の試料に放射線を照射して撮影することにより得られる減弱率と、前記試料の線減弱係数と厚さ及び前記放射線のスペクトルを用いて計算された減弱率とが一致するように、前記構成による放射線の吸収を計算するためのパラメータを補正することを特徴とする請求項１５に記載の画像処理装置。
17. 　前記構成による放射線の吸収を計算するための前記パラメータは、線源弱係数、厚み、充填率の少なくとも１つを含む、ことを特徴とする請求項１６に記載の画像処理装置。
18. 　補正対象のパラメータをユーザに選択させる選択手段をさらに備えることを特徴とする請求項１６または１７に記載の画像処理装置。
19. 　前記構成は、放射線エネルギーを可視光に変換する蛍光体を含むことを特徴とする請求項１５乃至１８のいずれか１項に記載の画像処理装置。
20. 　前記構成は、放射線発生装置と放射線撮像装置の間に配置される付加フィルタおよび／またはグリッドを含むことを特徴とする請求項１５乃至１９のいずれか１項に記載の画像処理装置。
21. 　前記補正手段は、前記エネルギーサブトラクション処理により算出された厚さが既知である試料の厚さが、前記試料の実際の厚さと一致するように前記放射線のスペクトルを補正することを特徴とする請求項１４に記載の画像処理装置。
22. 　前記補正手段は、前記算出された厚さと前記実際の厚さが一致するように、放射線発生装置と放射線撮像装置の間に配置された仮想の物体の減弱係数と厚みを計算し、前記仮想の物体の減弱係数と厚みを補正することにより前記放射線のスペクトルを補正することを特徴とする請求項２１に記載の画像処理装置。
23. 　前記補正手段は、前記計算された減弱率と前記実測により得られた減弱率が一致するように、放射線撮像装置のセンサの変換効率を計算し、前記変換効率を用いて前記放射線スペクトルを補正することを特徴とする請求項１４に記載の画像処理装置。
24. 　被写体に放射線を照射して撮影を行うことで得られた互いに異なる複数の放射線エネルギーに対応した複数の放射線画像を取得する取得工程と、
    　前記取得工程により取得された前記複数の放射線画像からエネルギーサブトラクション処理で用いられる複数の減弱率画像を生成する生成工程と、
    　前記放射線の線量、前記被写体の厚み、および前記放射線のエネルギーの少なくともいずれかに依存して生じる減弱率の誤差を低減するように、放射線画像または減弱率画像を補正する補正工程と、
    　前記補正工程による補正後の前記複数の減弱率画像を用いて前記エネルギーサブトラクション処理により物質特性画像を生成する処理工程と、を備えることを特徴とする画像処理方法。
25. 　エネルギーサブトラクション処理で用いられる放射線スペクトルを、計算により得られた画素値と実測により得られた画素値との差に基づいて補正する補正工程と、
    　被写体に放射線を照射して撮影を行うことで得られた互いに異なる複数の放射線エネルギーに対応した複数の放射線画像を取得する取得工程と、
    　前記取得工程により取得された複数の放射線画像から前記エネルギーサブトラクション処理に用いられる複数の減弱率画像を生成する生成工程と、
    　前記複数の減弱率画像と前記補正工程により補正された放射線スペクトルを用いて前記エネルギーサブトラクション処理により物質特性画像を生成する処理工程と、を備えることを特徴とする画像処理方法。
26. 　コンピュータを、請求項１乃至２３のいずれか１項に記載された画像処理装置の各手段として機能させるためのプログラム。

Images