(57)【要約】      （修正有）【課題】試料のＸ線吸収差により、画像再構成に生じる放射状パターンのアーチファクトを低減する。【解決手段】Ｘ線連続スペクトルを記憶し、単一スペクトルのＸ線投影データgを生成し、補正量を初期化し、Ｘ線投影データを補正量h、tで補正して撮影物に含む２つの物質のＸ線吸収係数のエネルギー依存性を考慮して撮影物に単一スペクトルのＸ線を照射した場合に相当する補正Ｘ線投影データg’を得、単一波長用画像再構成アルゴリズムを用いて逆投影演算してＣＴ画像を再構成し、各物質を区別する１つのしきい値TH1、TH2、…を用いてＣＴ画像から各物質の領域を抽出しこの領域のＸ線経路の長さを演算し、発生Ｘ線連続スペクトルのエネルギー、各物質の領域の長さ及び各物質のＸ線吸収係数を用いて連続スペクトルのＸ線投影データfを演算し、単一ペクトルのＸ線投影データgと連続スペクトルのＸ線投影データfとの差で補正量h、tを演算する。【選択図】図１
2011203160.tif 000002
【特許請求の範囲】
【請求項１】
発生Ｘ線連続エネルギースペクトルを既知情報として記憶する記憶段階と、
撮影対象物の単一エネルギースペクトルのＸ線投影データを生成するＸ線投影データ生成段階と、
補正量を初期化する初期化段階と、
前記単一エネルギースペクトルのＸ線投影データを前記補正量で補正して前記撮影対象物に含まれる少なくとも２つの物質のＸ線吸収係数のエネルギー依存性を考慮して前記撮影対象物に単一エネルギースペクトルのＸ線を照射した場合に相当する補正Ｘ線投影データを得る補正段階と、
前記補正Ｘ線投影データを単一波長用画像再構成アルゴリズムを用いて逆投影演算してＣＴ画像を再構成するＣＴ画像再構成段階と、
前記各物質を区別するための少なくとも１つのしきい値を用いて前記ＣＴ画像から前記各物質の領域を抽出して該領域のＸ線経路の長さを演算する長さ演算段階と、
前記発生Ｘ線連続エネルギースペクトルのエネルギー、前記各物質の領域の長さ及び前記各物質のＸ線吸収係数を用いて連続エネルギースペクトルのＸ線投影データを演算するＸ線投影データ演算段階と、
前記単一エネルギースペクトルのＸ線投影データと前記連続エネルギースペクトルのＸ線投影データとの差により前記補正量を演算する補正量演算段階と
を具備するＸ線ＣＴ画像再構成方法。
【請求項２】
前記物質が既知であり、さらに、前記補正段階、前記ＣＴ画像再構成段階、前記長さ演算段階、前記Ｘ線投影データ演算段階及び前記補正量演算段階を繰返して前記補正量が所定範囲に収束するようにする段階を具備する請求項１に記載のＸ線ＣＴ画像再構成方法。
【請求項３】
前記物質が未知であり、さらに、前記ＣＴ画像のヒストグラムから前記しきい値を推定するしきい値推定段階を具備する請求項１に記載のＸ線ＣＴ画像再構成方法。
【請求項４】
さらに、前記補正段階、前記ＣＴ画像再構成段階、前記しきい値推定段階、前記長さ演算段階、前記Ｘ線投影データ演算段階及び前記補正量演算段階を繰返して前記補正量が所定範囲に収束するようにする段階を具備する請求項３に記載のＸ線ＣＴ画像再構成方法。
【請求項５】
前記各物質のＸ線吸収係数は予め設定された請求項１に記載のＸ線ＣＴ画像再構成方法。
【請求項６】
さらに、
前記各物質kの係数akを
【数１】
2011203160.tif 000013
を用いて演算する段階と、
該演算された係数ak及び前記発生Ｘ線連続エネルギースペクトルのエネルギーEiを用いて前記物質のＸ線吸収係数μikを
【数２】
2011203160.tif 000014
により演算する段階と
を具備する請求項１に記載のＸ線ＣＴ画像再構成方法。
【請求項７】
さらに、
前記各物質kの係数akを
ak=αk・THn(n=1,2,…)+β
但し、THnは前記しきい値、αk、βkは定数
により演算する段階と、
該演算された係数ak及び前記発生Ｘ線連続エネルギースペクトルのエネルギーEiを用いて前記各物質kのＸ線吸収係数μikを
【数２】
2011203160.tif 000015
により演算する段階と
を具備する請求項１に記載のＸ線ＣＴ画像再構成方法。
【請求項８】
さらに、前記係数akにより前記物質kを推定する段階を具備する請求項６あるいは７に記載のＸ線ＣＴ画像再構成方法。
【請求項９】
発生Ｘ線連続エネルギースペクトルを既知情報として記憶する記憶手順と、
撮影対象物の単一エネルギースペクトルのＸ線投影データを生成するＸ線投影データ生成手順と、
補正量を初期化する初期化手順と、
前記単一エネルギースペクトルのＸ線投影データを前記補正量で補正して前記撮影対象物に含まれる少なくとも２つの物質のＸ線吸収係数のエネルギー依存性を考慮して前記撮影対象物に単一エネルギースペクトルのＸ線を照射した場合に相当する補正Ｘ線投影データを得る補正手順と、
前記補正Ｘ線投影データを単一波長用画像再構成アルゴリズムを用いて逆投影演算してＣＴ画像を再構成するＣＴ画像再構成手順と、
前記各物質を区別するための少なくとも１つのしきい値を用いて前記ＣＴ画像から前記各物質の領域を抽出して該領域のＸ線経路の長さを演算する長さ演算手順と、
前記発生Ｘ線連続エネルギースペクトルのエネルギー、前記各物質の領域の長さ及び前記各物質のＸ線吸収係数を用いて連続エネルギースペクトルのＸ線投影データを演算するＸ線投影データ演算手順と、
前記単一エネルギースペクトルのＸ線投影データと前記連続エネルギースペクトルのＸ線投影データとの差により前記補正量を演算する補正量演算手順と
を具備するＸ線ＣＴ画像再構成プログラム。
【請求項１０】
前記物質が既知であり、さらに、前記補正手順、前記ＣＴ画像再構成手順、前記長さ演算手順、前記Ｘ線投影データ演算手順及び前記補正量演算手順を繰返して前記補正量が所定範囲に収束するようにする手順を具備する請求項９に記載のＸ線ＣＴ画像再構成プログラム。
【請求項１１】
前記物質が未知であり、さらに、前記ＣＴ画像のヒストグラムから前記しきい値を推定するしきい値推定手順を具備する請求項９に記載のＸ線ＣＴ画像再構成プログラム。
【請求項１２】
さらに、前記補正手順、前記ＣＴ画像再構成手順、前記しきい値推定手順、前記長さ演算手順、前記Ｘ線投影データ演算手順及び前記補正量演算手順を繰返して前記補正量が所定範囲に収束するようにする手順を具備する請求項１１に記載のＸ線ＣＴ画像再構成プログラム。
【請求項１３】
前記各物質のＸ線吸収係数は予め設定された請求項９に記載のＸ線ＣＴ画像再構成プログラム。
【請求項１４】
さらに、
前記物質kの係数akを
【数１】
2011203160.tif 000016
を用いて演算する手順と、
該演算された係数ak及び前記発生Ｘ線連続エネルギースペクトルのエネルギーEiを用いて前記物質kのＸ線吸収係数μikを
【数２】
2011203160.tif 000017
により演算する手順と
を具備する請求項９に記載のＸ線ＣＴ画像再構成プログラム。
【請求項１５】
さらに、
前記物質kの係数akを
ak=αk・THn(n=1,2,…)+βk
但し、THnは前記しきい値、αk、βkは定数
により演算する手順と、
該演算された係数ak及び前記発生Ｘ線連続エネルギースペクトルのエネルギーEiを用いて前記各物質kのＸ線吸収係数μikを
【数２】
2011203160.tif 000018
により演算する手順と
を具備する請求項７に記載のＸ線ＣＴ画像再構成プログラム。
【請求項１６】
さらに、前記係数akにより前記物質kを推定する段階を具備する請求項１４あるいは１５に記載のＸ線ＣＴ画像再構成プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明はＸ線ＣＴ（Computed Tomography）画像再構成方法及びＸ線ＣＴ画像再構成プログラムに関する。
【背景技術】
【０００２】
一般に、Ｘ線ＣＴ画像再構成方法においては、撮影対象物の周囲を多数の方向から通常の波長分布を有するＸ線を投影し、Ｘ線の透過エネルギーとしての投影データを取得し、この投影データをフィルタ補正逆投影法、逐次近似法等の単一波長用画像再構成アルゴリズムを用いて逆投影して画像（ＣＴ画像）を再構成する。
【０００３】
上述のＸ線ＣＴ画像再構成方法において、撮影対象物にＸ線吸収係数が極端に大きい物質領域たとえば金属が存在すると、その物質領域を波長分布を有する照射Ｘ線が透過することによりその線質が変化するビームハードニング現象により、画像アーチファクトと呼ばれる放射状ノイズパターンがＣＴ画像に発生する。この画像アーチファクトは医療現場、非破壊検査現場等で再構成画像の判別の妨げとなる。
【０００４】
次に、ビームハードニング現象について、図４０、図４１、図４２、図４３を参照して説明する。
【０００５】
図４０は一般的なＸ線ＣＴ画像再構成装置を示す図である。
【０００６】
図４０において、１はＸ線発生器、２はＸ線検出器、３は撮影対象物、４は、Ｘ線発生器１、Ｘ線検出器２、撮影対象物３を制御する制御回路である。制御回路４はコンピュータで構成され、ＣＰＵ、ＲＯＭ、ＲＡＭ等を含む。
【０００７】
図４１は図４０の撮影対象物３の物質のＸ線吸収係数μを示す。
【０００８】
図４１に示すように、Ｘ線吸収係数μの分布が物質毎に異なるだけでなく、同一物質でもＸ線波長によってＸ線吸収係数μが異なる、つまり、エネルギー依存性を有する。
【０００９】
図４２は図４０のＣＴ画像再構成装置のＸ線エネルギースペクトルを示し、（Ａ)はＸ線発生器１の発生Ｘ線の連続エネルギースペクトル（以下、単に、連続スペクトル）を示し、（Ｂ）は大きさＬ、Ｘ線吸収係数μ（Ｅ）の撮影対象物３を透過後のＸ線検出器２の検出Ｘ線の連続スペクトルを示す。
【００１０】
図４２に示すように、撮影対象物３の透過前後でＸ線の連続スペクトルは異なる。
【００１１】
ところで、図４０のＸ線ＣＴ画像再構成装置におけるＣＴ画像再構成演算においては、Ｘ線吸収量のエネルギー依存性を考慮していない。つまり、μ一定の単一エネルギースペクトル（以下単に、単一スペクトル）のＸ線エネルギー照射つまり単一波長照射を前提としている。この理由は、Ｘ線連続スペクトルは装置、照射するＸ線等の種類、撮影条件等に依存するので、これに対応することが困難であるためである。しかし、実際のＣＴ撮影は通常のＸ線のために波長分布を有する。この結果、単一波長用ＣＴ画像再構成時の演算によって矛盾が発生し、図４３に示すような画像アーチファクトが発生する。
【００１２】
画像アーチファクトを低減する第１の従来のＸ線ＣＴ画像再生方法は投影データのシノグラムから高Ｘ線吸収係数の物質領域たとえば金属領域を抽出して仮想投影によりその物質領域の投影データを改変して修正投影データを演算し、この修正投影データを逆投影して再構成画像（ＣＴ画像）を再構成している（参照：特許文献１）。
【００１３】
上述の第１の従来のＸ線ＣＴ画像構成方法においては、厳しい画像アーチファクトを低減することができるが、軟組織部分の抽出が難しく投影データの修正が不充分であった。
【００１４】
画像アーチファクトを低減する第2の従来のＸ線ＣＴ画像構成方法は、投影条件を変更して複数回ＣＴ撮影を行って複数の投影データを得てＸ線連続スペクトルの種類を限定し、Ｘ線吸収量を求めている（参照：特許文献２、非特許文献１）。
【先行技術文献】
【特許文献】
【００１５】
【特許文献１】特開２００６−１６７１６１号公報
【特許文献２】特開２００８−２４１３７６号公報
【非特許文献】
【００１６】
【非特許文献１】李根旭、小関道彦、木村仁、伊能教夫：Ｘ線ＣＴ画像の画質向上に関する研究（Ｘ線スペクトルを考慮した画像再構成手法の提案）、計測自動制御学会  第24回センシングフォーラム資料、pp.107-111,2007
【発明の概要】
【発明が解決しようとする課題】
【００１７】
しかしながら、上述の第２の従来のＸ線ＣＴ画像再構成方法においては、複数のＸ線連続スペクトルを必要とするので、ＣＴ撮影回数が増加すると共に、撮影条件を変更することによる異なるＸ線連続スペクトルを得ることが難しいという課題がある。
【課題を解決するための手段】
【００１８】
上述の課題を解決するために本発明は、Ｘ線ＣＴ画像再構成方法において、発生Ｘ線連続エネルギースペクトルを既知情報として記憶し、撮影対象物の単一エネルギースペクトルのＸ線投影データを生成し、補正量を初期化し、単一エネルギースペクトルのＸ線投影データを補正量で補正して撮影対象物に含まれる少なくとも２つの物質のＸ線吸収係数のエネルギー依存性を考慮して撮影対象物に単一エネルギースペクトルのＸ線を照射した場合に相当する補正Ｘ線投影データを得、補正Ｘ線投影データを単一波長用画像再構成アルゴリズムを用いて逆投影演算してＣＴ画像を再構成し、各物質を区別するための少なくとも１つのしきい値を用いてＣＴ画像から各物質の領域を抽出してこれらの領域のＸ線経路の長さを演算し、発生Ｘ線連続エネルギースペクトルのエネルギー、各物質の領域の長さ及び各物質のＸ線吸収係数を用いて連続エネルギースペクトルのＸ線投影データを演算し、単一エネルギースペクトルのＸ線投影データと連続エネルギースペクトルのＸ線投影データとの差により補正量を演算するものである。すなわち、従来の単一波長用の画像再構成アルゴリズムを用いて得られるＣＴ画像（元ＣＴ画像）に対して、推定される撮影対象物に対して単一エネルギースペクトルのＸ線を照射したときＸ線投影データを推定する。この推定には、元ＣＴ画像に少なくとも１つのしきい値を用いて少なくとも２つの物質の領域を抽出し、各Ｘ線経路の長さを演算する。次いで、既知情報のＸ線連続エネルギースペクトルのエネルギー、物質の領域の長さ及びＸ線吸収係数を用いて実際に照射された連続エネルギースペクトルのＸ線投影データを補正演算し、単一エネルギースペクトルのＸ線投影データに近い値を得るようにしたものである。このようにして、少なくとも１つのしきい値で区分された少なくとも２つの物質が及ぼすビームハードニングの度合を個々に換算して補正量で表わす。
【００１９】
また、上述の補正から補正量演算を繰返して補正量が所定範囲内に収束させる。この場合、しきい値が未知であれば、ＣＴ画像のヒストグラムによりしきい値を推定してしきい値の数を増加させる。このようにして、補正量が更新され、画像アーチファクトが徐々に減少することで物質の領域のＸ線経路の長さが正確となる。
【発明の効果】
【００２０】
本発明によれば、単一エネルギースペクトルのＸ線投影データは変更されないので、ＣＴ撮影は１回で済み、しかも、補正を繰返すことによって画像アーチファクトを低減することができる。
【図面の簡単な説明】
【００２１】
【図１】本発明に係るＸ線ＣＴ画像再構成方法の第１の実施の形態を説明するためのフローチャートである。
【図２】図１のステップ１０１のＸ線連続スペクトルの例を示すグラフである。
【図３】図１のステップ１０６のしきい値を説明する図であって、（A）はＣＴ画像、（B）は画素値ヒストグラムである。
【図４】図１のステップ１０８の補正量を説明するためのグラフである。
【図５】Ｘ線発生器の発生Ｘ線の連続スペクトルを示す表である。
【図６】生体骨のＸ線吸収係数を示す表である。
【図７】図１の第１のシミュレーション例に用いられる撮影対象物としての骨円柱モデルを示す図である。
【図８】図７の補正前の骨円柱モデルを示し、（A）はＣＴ画像、（B）は（A）センタライン上のＣＴ値を示すグラフである。
【図９】図８の（A）のＣＴ画像から物質の領域抽出を示すＣＴ画像である。
【図１０】図７の1回目補正の骨円柱モデルを示し、（A）はＣＴ画像、（B）は（A）のセンタライン上のＣＴ値を示すグラフである。
【図１１】鉄のＸ線吸収係数を示す表である。
【図１２】アクリル樹脂のＸ線吸収係数を示す表である。
【図１３】図１の第２のシミュレーション例に用いられる撮影対象物としての複数同種ピンモデルを示す図であって、（A)は斜視像、（B）は上面像である。
【図１４】図１３の補正前の複数同種ピンモデルを示し、（A）はＣＴ画像、（B）は（A）センタライン上のＣＴ値を示すグラフである。
【図１５】図１４の（A）のＣＴ画像から物質の領域抽出を示し、（A）は鉄の領域抽出を示すＣＴ画像、（B）はアクリル樹脂の領域抽出を示すＣＴ画像である。
【図１６】図１３の1回目補正の複数同種ピンモデルを示し、（A）はＣＴ画像、（B）は（A）のセンタライン上のＣＴ値を示すグラフである。
【図１７】アルミニウムのＸ線吸収係数を示す表である。
【図１８】図１の第３のシミュレーション例に用いられる撮影対象物としての複数異種ピンモデルを示す上面像である。
【図１９】図１８の補正前の複数異種ピンモデルを示し、（A）はＣＴ画像、（B）は（A）センタライン上のＣＴ値を示すグラフである。
【図２０】図１９の画素値ヒストグラムを示すグラフである。
【図２１】図１９の（A）のＣＴ画像から物質の領域抽出を示し、（A）は鉄の領域抽出を示すＣＴ画像、（B）はアルミニウムの領域抽出を示すＣＴ画像、（C）はアクリル樹脂の領域抽出を示すＣＴ画像である。
【図２２】図１８の１回目補正の複数異種ピンモデルを示し、（A）はＣＴ画像、（B）は（A）の画素値ヒストグラムを示すグラフである。
【図２３】図１８の２回目補正の複数異種ピンモデルを示し、（A）はＣＴ画像、（B）は（A）の画素値ヒストグラムを示すグラフである。
【図２４】図１８の３回目補正の複数異種ピンモデルを示し、（A）はＣＴ画像、（B）は（A）の画素値ヒストグラムを示すグラフである。
【図２５】図１の第３のシミュレーションの補正量の収束を示すグラフである。
【図２６】図１の第１の実験例を示し、（A）は補正前のＣＴ画像、（B）は補正後のＣＴ画像である。
【図２７】図１の第２の実験例を示し、（A）は補正前のＣＴ画像、（B）は補正後のＣＴ画像である。
【図２８】図１の第３の実験例を示し、（A）は鉄釘―米粒モデルの斜視像、（B)は撮影時のＸ線スペクトルを示すグラフである。
【図２９】図１の第３の実験例を示し、（A）は補正前のＣＴ画像、（B）は補正後のＣＴ画像である。
【図３０】第３の実験例において鉄のみ抽出による補正後のＣＴ画像である。
【図３１】本発明に係るＸ線ＣＴ画像再構成方法の第２の実施の形態を説明するためのフローチャートである。
【図３２】図３１のステップ２０２、２０３における物質鉄、アルミニウム、アクリル樹脂のＸ線吸収係数の例を示すグラフである。
【図３３】図３１の第１のシミュレーション例の1回目補正後のＣＴ画像である。
【図３４】図３１の第２のシミュレーション例を示し、（A）は1回目補正後のＣＴ画像、（B）は補正後のＣＴ画像である。
【図３５】本発明に係るＸ線ＣＴ画像再構成方法の第３の実施の形態を説明するためのフローチャートである。
【図３６】図３５の第１のシミュレーション例を示し、（A）、（B）、（C）、（D）、（E）、（F）は、補正前、1回目補正後、2回目補正後、3回目補正後、4回目補正後、5回目補正後のＣＴ画像である。
【図３７】図３５の第２のシミュレーション例を示し、（A）、（B）、（C）、（D）、（E）は、補正前、1回目補正後、2回目補正後、3回目補正後、4回目補正後のＣＴ画像である。
【図３８】図３５の第１のシミュレーション例に用いられる撮影対象物としての生体組織モデルを示す図である。
【図３９】図３８のＣＴ画像を示し、（A）は補正前のＣＴ画像、（B）は補正後のＣＴ画像である。
【図４０】一般的なＸ線ＣＴ画像再構成装置を示す図である。
【図４１】図４０の撮影対象物の物質のＸ線吸収係数を示すグラフである。
【図４２】図４０のＸ線連続スペクトルを示すグラフである。
【図４３】画像アーチファクトの例を示すＣＴ画像である。
【発明を実施するための形態】
【００２２】
図１は本発明に係るＸ線ＣＴ画像再構成方法の第１の実施の形態を説明するためのフローチャートである。図１のフローチャートは図４０のＸ線ＣＴ画像再構成装置によって実行される。
【００２３】
始めに、ステップ１０１では、Ｘ線発生器１のＸ線連続スペクトルを予め記憶する。このＸ線発生器１のＸ線管の種類、撮影条件によって定まるＸ線連続スペクトルは専用の機器で予め計測できる。尚、Ｘ線連続スペクトルは実際には離散値としてEi及びその光子数ｎ（Ei）で記憶される。
【００２４】
次に、ステップ１０２では、後述の補正量ht及び補正回数tを初期化する。つまり、
h0←0
t←0
ここで、h0は投影データの数だけ存在するので、h0はその総称を表している。
【００２５】
次にステップ１０３では、Ｘ線検出器２の検出信号により単一スペクトルつまり波長iでの単一スペクトル投影データ（Ｘ線吸収量）gを生成する。尚、波長iはＸ線管電圧から計算される代表値を採用する。
【００２６】
ステップ１０３におけるＸ線吸収量gを理論的に考察すると、Ｘ線吸収量gは数１の式で表される。
【数１】
2011203160.tif 000003
但し、Iin:物質kの入射時のエネルギー
Iout：物質kの透過後のエネルギー
Lk：物質kのＸ線経路の長さ
μik：波長i及び物質kに依存するＸ線吸収係数
Ｘ線が透過する物質がn種類の物質で構成されていれば、数１の式は数２の式となる
【数２】
2011203160.tif 000004
数２の式を対数で表すと、数３の式となる。
【数３】
2011203160.tif 000005
【００２７】
次に、ステップ１０４では、単一スペクトル投影データgに補正量htを減算して補正投影データg’を求める。
g’= g - ht
尚、ここでも、g’は投影データの数だけ存在するので、その総称を表している。
【００２８】
次に、ステップ１０５では、補正投影データg’をフィルタ補正逆投影法等により逆投影してＣＴ画像を再構成する。この場合の逆投影法は、従来の単一スペクトルのＸ線投影データを再構成するための単一波長用画像再構成アルゴリズムを用いる。
【００２９】
次に、ステップ１０６では、既知のしきい値TH1、TH2、…を用いてＣＴ画像から物質kのＸ線経路の領域の長さLkを演算する。たとえば、図３の（Ａ）に示すＣＴ画像が図３の（Ｂ）に示すごとく３種類の物質、つまり、鉄、アルミニウム、アクリル樹脂を含んでいれば、ステップ１０６において、TH1を画素値14000、TH2を画素値600に予め設定しておく。尚、対象撮影物３が複数の未知の物質を含んでいる場合には、ステップ１０６において、図３の（Ｂ）の画素値ヒストグラムを演算し、画素値の集中箇所、たとえば、図３の（Ｂ）における画素値-900、1800、8000、20000の集中箇所を判別し、その間にしきい値を推定して設定する。尚、このしきい値推定アルゴリズムは容易に得ることができる。
【００３０】
次に、ステップ１０７では、各単一スペクトル投影データgに対応する連続スペクトル投影データfを演算する。
【００３１】
連続スペクトル投影データfはステップ１０１にて記憶されたＸ線発生器１のＸ線連続スペクトルのエネルギーEi、物質kのＸ線吸収係数μik及びステップ１０６にて演算された物質kのＸ線経路の長さLkに基づいて演算される。つまり、物質kの入射時のＸ線エネルギーIinは数４の式で表わされる。
【数４】
2011203160.tif 000006
また、物質kの透過後のＸ線エネルギーIoutは数５の式で表わされる。
【数５】
2011203160.tif 000007
数４、数５の式から数６の式が得られる。
【数６】
2011203160.tif 000008
数６の式を対数で表わすと、連続スペクトル投影データ（Ｘ線吸収量）fは数７の式で表わされる。
【数７】
2011203160.tif 000009
【００３２】
次に、ステップ１０８では、
t ← t + 1
とする。
【００３３】
次に、ステップ１０９では、t回目の補正量htを演算する。すなわち、
ht = g - f
【００３４】
すなわち、図４に示すように、物質kのＸ線経路の長さLkが同一であっても、単一スペクトル投影データ（Ｘ線吸収量）gと連続スペクトル投影データ（Ｘ線吸収量）fとでは値が補正量htだけ異なっている。この補正量htだけ異なっている分、この補正量htが画像アーチファクトとして現われている。ステップ１０６における画像アーチファクトが存在する投影データg’からの物質kのＸ線経路の長さLkにも僅差が存在する。従って、ステップ１０４において単一スペクトル投影データgをこの補正量htで減算して補正することにより連続スペクトル投影データfに近づけることができ、長さLkも安定的収束し、これに伴い、画像アーチファクト低減効果が生ずる。
【００３５】
ステップ１１０はステップ１０４〜１０９のフローを補正量htが所定範囲に収束して|ht-ht-1|がε未満になるまで繰返す。|ht-ht-1|がε未満となったときにステップ１１１にてこのルーチンは終了する。尚、ここでも、|ht-ht-1|は投影データの数だけ存在し、その総和を意味する。
【００３６】
図１のフローチャートで重要なのは、1回のＣＴ撮影で得られた単一スペクトル投影データgは変更されず、その補正量htのみを更新し、この更新によって画像アーチファクトを徐々に低減することで物質kのＸ線経路の長さLkの演算も正確になっていくことである。
【００３７】
図１のフローチャートによるシミュレーション例を説明する。つまり、この場合、制御回路（コンピュータ）４上でＣＴ撮影する仮想投影手法を用い、Ｘ線発生器１のＸ線連続スペクトルは図５に示す９個の離散値で近似する。
【００３８】
第１のシミュレーション例は図６に示すＸ線吸収係数μikを有する物質つまり生体骨を用いた図７に示す骨円柱モデルに適用する。つまり、物質k及びそのＸ線吸収係数μikは既知とする。この場合、骨円柱モデルは均一であり、かつ形状が単純であるので、ビームハードニングの影響が小さく、従って、画像アーチファクトが少ないと考えられるケースである。
【００３９】
第１のシミュレーション例においては、ステップ１０１において、図５のＸ線連続スペクトル及び図６のＸ線吸収係数μikを記憶する。そして、ステップ１０３において、これらに基づいて単一スペクトル投影データgを生成し、ステップ１０５において、単一スペクトル投影データgからＣＴ画像を再構成する。この場合、図８の（Ａ）に示すＣＴ画像が得られ、一見すると画像アーチファクトは出現していないが、図８の（Ｂ）に示すように、ＣＴ画像のセンタライン上のＣＴ値をプロットすると、円柱中央部の画素で輝度が上昇し、円柱周辺部の画素で輝度が低下していることが分かる。つまり、円柱周辺部の画素のＣＴ値が真値1240より大きくなっている。このように、単純な形状でもビームハードニングの影響で適切なＸ線吸収量が得られず、ＣＴ値の定量性が喪失してしまう。
【００４０】
次に、図８の（Ａ）のＣＴ画像に１回目補正を行う。つまり、ステップ１０６において、図８の（Ａ）のＣＴ画像の画素に対して既知のしきい値TH1たとえば500を用いて図９に示す物質生体骨の領域を抽出する。図９に示すごとく、物質生体骨の領域は正確に抽出されている。この抽出領域のＸ線経路毎に長さLkを演算する。次いで、ステップ１０７において、その長さLk、図５のＸ線連続スペクトルのエネルギーEi、図６のＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、１回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データから図１０の（Ａ）に示すＣＴ画像を再構成する。
【００４１】
図１０の（Ａ）の１回目補正のＣＴ画像についてセンタライン上の画素をプロットすると、図１０の（Ｂ）に示すごとく、円柱周辺部の画素の輝度が上昇して円柱内部で均一となり、50keVのＣＴ真値1240に近づいていることが分かる。
【００４２】
第２のシミュレーション例は図１１、図１２に示すＸ線吸収係数μikを有する物質つまり鉄、アクリル樹脂を用いた図１３に示す複数同種ピンモデルに適用する。類似モデルとしては、電子部品、半導体の非破壊検査の際に、金属部品が樹脂中に埋め込められていると、金属部品の周囲に画像アーチファクトが発生して金属部品の周囲の樹脂部分の形状、組成が判別できなくなる場合がある。また、人体にインプラント金属が埋め込められていると、インプラント金属の周囲に画像アーチファクトが発生してインプラント金属の周囲の筋肉部分の形状、組成が判別できなくなる場合がある。尚、この場合も、物質k及びそのＸ線吸収係数μikは既知とする。この複数同種ピンモデルは形状が少し複雑であるので、ビームハードニングの影響が大きく、従って、画像アーチファクトがある程度大きいと考えられるケースである。
【００４３】
第２のシミュレーション例においては、ステップ１０１において、図５のＸ線連続スペクトルのエネルギーEi及び図１１、図１２のＸ線吸収係数μikを記憶する。そして、ステップ１０３において、これらに基づいて単一スペクトル投影データgを生成し、ステップ１０５において、単一スペクトル投影データgからＣＴ画像を再構成する。この場合、図１４の（Ａ）に示すＣＴ画像が得られ、強い画像アーチファクトが出現している。図１４の（Ｂ）に示すように、ＣＴ画像のセンタライン上のＣＴ値をプロットすると、ピン中央部の画素で輝度が急上昇し、ピン周辺部の画素で輝度が急低下していることが分かる。つまり、ピン周辺部の画素のＣＴ値が真値20000より大きくなっている。このように、複雑な形状では強いビームハードニングの影響で適切なＸ線吸収量が得られず、ＣＴ値の定量性が喪失してしまう。
【００４４】
次に、図１４の（Ａ）のＣＴ画像に１回目補正を行う。つまり、ステップ１０６において、図１４の（Ａ）のＣＴ画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて図１５に示す物質鉄及びアクリル樹脂の各領域を抽出する。この場合、図１５の（Ａ）に示すごとく、物質鉄の領域は正確に抽出されているが、図１５の（Ｂ）に示すごとく、物質アクリル樹脂の領域は物質鉄に起因するアーチファクトの影響で不正確である。これらの抽出領域のＸ線経路毎に長さLkを演算する。次いで、ステップ１０７において、これらの長さLk、図５のＸ線連続スペクトルのエネルギーEi、図１１、図１２のＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、１回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図１６の（Ａ）に示すＣＴ画像を再構成する。
【００４５】
図１６の（Ａ）の１回目補正のＣＴ画像についてセンタライン上画素をプロットすると、図１６の（Ｂ）に示すごとく、ピン周辺部の画素の輝度が上昇してピン内部で均一となり、50keVのＣＴ真値20000に近づいていることが分かる。また、同時に、物質アクリル樹脂の領域の画素の輝度も均一となり、50keVのＣＴ真値-680に近づいていることが分かる。
【００４６】
第３のシミュレーション例は図１１、図１２、図１７に示すＸ線吸収係数μikを有する物質つまり鉄、アクリル樹脂、アルミニウムを用いた図１８に示す複数異種ピンモデルに適用する。尚、この場合も、物質k及びそのＸ線吸収係数μikは既知とする。この複数異種ピンモデルは複数の高Ｘ線吸収係数を有する物質を有するので、ビームハードニングの影響が非常に大きく、従って、画像アーチファクトが大きいと考えられるケースである。
【００４７】
第３のシミュレーション例においては、ステップ１０１において、図５のＸ線連続スペクトルのエネルギーEi及び図１１、図１２、図１７のＸ線吸収係数μikを記憶する。そして、ステップ１０３において、これらに基づいて単一スペクトル投影データgを生成し、ステップ１０５において、単一スペクトル投影データgからＣＴ画像を再構成する。この場合、図１９の（Ａ）に示すＣＴ画像が得られ、物質鉄に起因する強い画像アーチファクトが出現している。図１９の（Ｂ）に示すように、ＣＴ画像のセンタライン上のＣＴ値をプロットすると、ピン中央部の画素で輝度が急上昇し、ピン周辺部の画素で輝度が急低下していることが分かる。つまり、ピン周辺部の画素のＣＴ値が真値20000あるいは1800より大きくなっている。このときの画素ヒストグラムを図２０に示すと、物質鉄の分布は真値20000から離れた23000程度の箇所にある。従って、図１のステップ１０６におけるしきい値TH1を8000と推定して設定し、物質鉄の領域を抽出できる。尚、この場合、物質鉄が未知であっても、１つのしきい値TH1が推定され設定される。この結果、図２１の（Ａ）に示すごとく、物質鉄の領域を正確に抽出できるが、図２１の（Ｂ）、（Ｃ）に示すごとく、物質アルミニウム、アクリル樹脂の各領域は物質鉄に起因するアーチファクトの影響で正確には抽出できない。
【００４８】
次に、図１９の（Ａ）のＣＴ画像に１回目補正を行う。つまり、上述の抽出領域のＸ線経路毎に長さLkを演算する。次いで、これらの長さLk、図５のＸ線連続スペクトルのエネルギーEi、図１１、図１２、図１７のＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、１回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図２２の（Ａ）に示すＣＴ画像を再構成する。図２２の（Ａ）のＣＴ画像においては、画像アーチファクトが残り、補正は不充分である。この理由は、物質鉄に起因する画像アーチファクトの画素値が物質アルミニウムの画素値とほぼ同一であるために、画像アーチファクトの形状が物質アルミニウムが存在するものとして再構成されたものと考えられる。このときの画素ヒストグラムを図２２の（Ｂ）に示すと、物質アルミニウム及びアクリル樹脂のピーク値がこれらの真値-680及び1800に現われている。従って、図１のステップ１０６におけるしきい値TH2を-680と1800との間たとえば1200と推定して設定すれば、物質アルミニウム、アクリル樹脂の各領域は正確に抽出できる。尚、この場合も、物質アルミニウム、アクリル樹脂が未知であっても、しきい値TH2が推定され設定される。従って、撮影対象物に含まれるｎ個の物質が未知である場合には、図１のステップ１０４〜１１０を繰返すことによりｎ−１個のしきい値が推定され設定されることになる。この結果、未知の物質の領域も正確に抽出できる。
【００４９】
次に、図２２の（Ａ）のＣＴ画像に２回目補正を行う。つまり、ステップ１０６において、図２２の（Ａ）のＣＴ画像の画素に対して推定のしきい値TH1たとえば8000及びしきい値TH2たとえば1200を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。次いで、ステップ１０７において、これらの長さLk、図５のＸ線連続スペクトルのエネルギーEi、図１１、図１２、図１７のＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、２回目の補正量h2を演算し、補正量h2で単一スペクトル投影データgを補正し、補正された投影データg’から図２３の（Ａ）に示すＣＴ画像を再構成する。
【００５０】
次に、図２３の（Ａ）のＣＴ画像に３回目補正を行う。つまり、ステップ１０６において、図２３の（Ａ）のＣＴ画像の画素に対して推定のしきい値TH1たとえば8000及びしきい値TH2たとえば1200を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。次いで、ステップ１０７において、これらの長さLk、図５のＸ線連続スペクトルのエネルギーEi、図１１、図１２、図１７のＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、３回目の補正量h3を演算し、補正量h3で単一スペクトル投影データgを補正し、補正された投影データg’から図２３の（Ａ）に示すＣＴ画像を再構成する。
【００５１】
図２４の（Ａ）の３回目補正のＣＴ画像についてセンタライン上画素をプロットすると、図２４の（Ｂ）に示すごとく、ピン周辺部の画素の輝度が上昇してピン内部で均一となり、50keVのＣＴ真値20000もしくは1800に近づいていることが分かる。また、同時に、物質アクリル樹脂の領域の画素の輝度も均一となり、50keVのＣＴ真値-680に近づいていることが分かる。
【００５２】
第３のシミュレーション例においては、ＣＴ画像をたとえ４回目補正、５回目補正、…を行ってもＣＴ画像はほとんど変わらない。つまり、補正回数が増加するにつれて、画像アーチファクトが低減し、物質鉄、アルミニウム、アクリル樹脂の領域がほぼ正確に抽出できる。従って、図２５に示すごとく、補正回数が３回を超えて総和|ht-ht-1|がε以下となると、各物質がその真値に近づいており、ＣＴ画像の誤差は低減し、補正の効果が向上していることが確認できる。このとき、補正量と前回補正量との差|ht-ht-1|も0に近づき、補正量が一定値に収束していることが分かる。
【００５３】
次に、図１のフローチャートによる実験例を説明する。
【００５４】
図２６は上述の第２のシミュレーション例に対応して行った第１の実験例である。撮影時の発生Ｘ線の連続スペクトルは不明であったので、図５の連続スペクトルの離散値で近似させた。また、物質鉄、アクリル樹脂のＸ線吸収係数μikは図１１、図１２に示す値を既知として用いた。さらに、物質鉄、アクリル樹脂の各領域の抽出のためのしきい値に適当な値を既知として用いた。この結果、補正前のＣＴ画像は図２６の（Ａ）に示すごとくであったが、図１のフローチャートに基づく補正を行った結果、図２６の（Ｂ）に示すＣＴ画像が得られた。つまり、画像アーチファクトは大幅に減少し、補正の効果が認められた。
【００５５】
図２７は鉄の十字架形状モデルに対して行った第２の実験例である。この場合も、撮影時の発生Ｘ線の連続スペクトルは不明であったので、図５の連続スペクトルの離散値で近似させた。また、物質鉄、アクリル樹脂のＸ線吸収係数μikは図１１、図１２に示す値を既知として用いた。さらに、物質鉄、アクリル樹脂の各領域の抽出のためのしきい値に適当な値を既知として用いた。この結果、補正前のＣＴ画像は図２７の（Ａ）に示すごとくであったが、図１のフローチャートに基づく補正を行った結果、図２７の（Ｂ）に示すＣＴ画像が得られた。やはり、画像アーチファクトは大幅に減少し、補正の効果が認められた。
【００５６】
図２８の（Ａ）に示す９本の鉄釘を米粒の山へ埋込んだモデルに対して行ったものである。撮影時の管電圧130kVの発生Ｘ線の連続スペクトルは、図２８の（Ｂ）の連続スペクトルの離散値で近似させた。また、物質鉄、米粒のＸ線吸収係数μikは図１１、図１２に示す値を既知として用いた。つまり、米粒はアクリル樹脂とみなした。さらに、物質鉄、米粒の各領域の抽出のためのしきい値に13500及び750を既知として用いた。この結果、補正前のＣＴ画像は図２９の（Ａ）に示すごとくであった。つまり、鉄釘と鉄釘とを結ぶラインで強い画像アーチファクトが発生している。これに対し、図１のフローチャートに基づく補正を行った結果、図２９の（Ｂ）に示すＣＴ画像が得られた。やはり、画像アーチファクトは大幅に減少し、補正の効果が認められた。
【００５７】
尚、第３の実験例において、しきい値13500だけ設定して鉄釘（鉄）のみの領域を抽出した場合には、図３０に示すごとく、米粒が正確に抽出できなかった。この理由は、連続スペクトルにより、鉄釘（鉄）だけでなく、米粒部分でもＸ線吸収量の過吸収が発生しているものと考えられるからである。
【００５８】
このように、実験例によれば、本発明はＣＴ画像を用いた個体別モデリング、リバースエンジニアリング非破壊検査で有効であることが分かる。
【００５９】
図３１は本発明に係るＸ線ＣＴ画像再構成方法の第２の実施の形態を説明するためのフローチャートであって、対象撮影物３を構成する物質kの数が少なくかつこれら物質kのＸ線吸収係数μikが未知の場合に適用される。
【００６０】
図３１においては、図１において、ステップ１０１をステップ２０１に置換し、ステップ２０２、２０３を追加したものである。すなわち、ステップ２０１では、図５のＸ線連続スペクトルのエネルギーEiのみを記憶し、ステップ２０２、２０３において、Ｘ線連続スペクトルのエネルギーEi及び物質kのＸ線経路の長さLkに基づいて物質kのＸ線吸収係数μikを演算する。
【００６１】
ステップ２０２、２０３の根拠を説明する。Ｘ線検出器２のある位置の投影データDは数８の式で表わされる。
【数８】
2011203160.tif 000010
ここで、Ｘ線連続スペクトルのエネルギーEiが既知であるので、数８の式の投影データはＸ線吸収係数μik及びステップ１０６で定まる長さLkの関数となる。他方、Ｘ線吸収係数μikは数９の式で表わされる。
【数９】
2011203160.tif 000011
数９の式の妥当性は、図３２に示すように、物質鉄、アルミニウム、アクリル樹脂に対してak=134000、7000、200としたときに、数９の式により求められたＸ線吸収係数μikは実際のＸ線吸収係数μikにほぼ一致することから確認できる。
数９の式で数８の式を表わすと、数１０の式となる。
【数１０】
2011203160.tif 000012
数１０の式において、物質kの種類が未知数ｎ個のｎ次方程式としてakを解くことができる。実際に検出されるデータはＸ線検出器２の列数たとえば512×撮影の回転方向数たとえば800個分入手できるので、物質kの種類数nはどんなに多くとも解くことができる。つまり、複数の方向から照射したＸ線ビーム毎に投影データD’が得られるので、数１０の式がその数だけ成立する。この場合、数１０の式の数はｎ以上であればakを求めることができる。但し、領域の長さLkは画像アーチファクトによる誤差を有するので、akは最小２乗法で求める。尚、このとき、求められたakに基づき、物質kを推定できる。
【００６２】
従って、ステップ２０２では、数１０の式のakを最小２乗法により演算する。次いで、ステップ２０３において、物質kのＸ線吸収係数μikを数９の式により演算する。
【００６３】
図３１のフローチャートによるシミュレーション例を説明する。つまり、この場合、制御回路（コンピュータ）４上でＣＴ撮影する仮想投影手法を用い、Ｘ線発生器１のＸ線連続スペクトルは図５に示す９個の離散値で近似する。
【００６４】
第１のシミュレーション例は２つの物質つまり鉄、アクリル樹脂を用いた図１３に示す複数同種ピンモデルに適用する。但し、物質鉄、アクリル樹脂のＸ線吸収係数μikは未知とする。
【００６５】
第１のシミュレーション例においては、ステップ２０１において、図５のＸ線連続スペクトルのエネルギーEiを記憶する。そして、ステップ１０３において、これらに基づいて単一スペクトル投影データgを生成し、ステップ１０５において、単一スペクトル投影データgからＣＴ画像を再構成する。
【００６６】
次に、ＣＴ画像に１回目補正を行う。つまり、ステップ１０６において、ＣＴ画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のＸ線経路毎に長さLkを演算する。次いで、ステップ２０２において、数１０の式のakを最小２乗法により演算する。その結果、
a1 = 124651
a2 = 41
が得られた。尚、a（鉄）、a（アクリル樹脂）の真値は130000、200であるので、上述の結果値は真値に近いことが分かる。従って、a1 = 124651は鉄に相当、a2= 41はアクリル樹脂に相当すると推定する。次いで、ステップ２０３において、物質鉄及びアクリル樹脂相当のＸ線吸収係数μikを数９の式を用いて演算する。次いで、ステップ１０７において、上述の長さLk、図５のＸ線連続スペクトルのエネルギーEi、ステップ２０３にて得られたＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、１回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図３３に示すＣＴ画像を再構成する。このように、Ｘ線吸収係数μikが未知である場合にも、画像アーチファクトが大幅に減少する。
【００６７】
第２のシミュレーション例は３つの物質つまり鉄、アルミニウム及びアクリル樹脂を用いた図１８に示す複数異種ピンモデルに適用する。但し、物質鉄、アルミニウム及びアクリル樹脂のＸ線吸収係数μikは未知とする。
【００６８】
第２のシミュレーション例においては、ステップ２０１において、図５のＸ線連続スペクトルのエネルギーEiを記憶する。そして、ステップ１０３において、これらに基づいて単一スペクトル投影データgを生成し、ステップ１０５において、単一スペクトル投影データgからＣＴ画像を再構成する。
【００６９】
次に、ＣＴ画像に１回目補正を行う。つまり、ステップ１０６において、ＣＴ画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。これらの抽出領域のＸ線経路毎に長さLkを演算する。次いで、ステップ２０２において、数１０の式のakを最小２乗法により演算する。その結果、
a1 = 127682
a2 = -1702
a3 = 416
が得られた。この場合、画像アーチファクトの影響によりアルミニウム相当の領域抽出が不適切であるために、a（アルミニウム）の真値7000から離れてしまった。また、Ｘ線吸収係数は負とならないので、
a1 = 127682
a2 = 0
a3 = 416
とする。次いで、ステップ２０３において、物質鉄、アルミニウム及びアクリル樹脂相当のＸ線吸収係数μikを数９の式を用いて演算する。次いで、ステップ１０７において、上述の長さLk、図５のＸ線連続スペクトルのエネルギーEi、ステップ２０３にて得られたＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、１回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図３４の（Ａ）に示すＣＴ画像を再構成する。このように、１回目補正により、不充分であるが、画像アーチファクトが減少する。
【００７０】
次に、図３４の（Ａ）のＣＴ画像に２回目補正を行う。つまり、ステップ１０６において、ＣＴ画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。これらの抽出領域のＸ線経路毎に長さLkを演算する。次いで、ステップ２０２において、数１０の式のakを最小２乗法により演算する。その結果、
a1 = 123476
a2 = 6962
a3 = 21
が得られた。尚、a（鉄）、a（アルミニウム）、a（アクリル樹脂）の真値は130000、7000、200であるので、上述の結果値は真値に近いことが分かる。従って、a1 = 123476は鉄、a2= 6962はアルミニウム、a3 = 21はアクリル樹脂に相当すると推定する。次いで、ステップ２０３において、物質鉄、アルミニウム及びアクリル樹脂相当のＸ線吸収係数μikを数９の式を用いて演算する。次いで、ステップ１０７において、上述の長さLk、図５のＸ線連続スペクトルのエネルギーEi、ステップ２０３にて得られたＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、１回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図３４の（Ｂ）に示すＣＴ画像を再構成する。このように、２回目補正により、画像アーチファクトが大幅に減少する。
【００７１】
図３５は本発明に係るＸ線ＣＴ画像再構成方法の第３の実施の形態を説明するためのフローチャートであって、対象撮影物３を構成する物質kの数が多くかつこれら物質kのＸ線吸収係数μikが未知の場合に適用される。
【００７２】
図３５においては、図３１において、ステップ２０２をステップ２０２Ａに置換したものである。すなわち、ステップ２０１Ａでは、図３１のフローチャートにおけるakを
ak = 9.873・THn + 9873
= 9.873・(THn+ 1000)
但し、THn(n=1,2,…)は暫定的なしきい値、
により演算する。ここで、値9.873、9873は、空気のＣＴ値が-1000である場合にa（空気）=0、鉄のＣＴ値12572の場合にa（鉄）=134000となるように線形に規定したものである。ここで、鉄のＣＴ値12572は実際のＣＴ画像の最大ＣＴ値である。但し、9.873、9873は与えられた2点が線形となるように、他の値αk、βkとなし得る。この場合も、求められたakに基づき、物質kを推定できる。
【００７３】
図３５のフローチャートによるシミュレーション例を説明する。つまり、この場合、制御回路（コンピュータ）４上でＣＴ撮影する仮想投影手法を用い、Ｘ線発生器１のＸ線連続スペクトルは図５に示す９個の離散値で近似する。
【００７４】
第１のシミュレーション例は２つの物質つまり鉄、アクリル樹脂を用いた図１３に示す複数同種ピンモデルに適用する。但し、物質鉄、アクリル樹脂のＸ線吸収係数μikは未知とする。
【００７５】
第１のシミュレーション例においては、ステップ２０１において、図５のＸ線連続スペクトルのエネルギーEiを記憶する。そして、ステップ１０３において、これらに基づいて単一スペクトル投影データgを生成し、ステップ１０５において、単一スペクトル投影データgから図３６の（Ａ）に示すＣＴ画像を再構成する。
【００７６】
次に、ＣＴ画像に１回目補正を行う。つまり、ステップ１０６において、ＣＴ画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のＸ線経路毎に長さLkを演算する。次いで、ステップ２０２Ａにおいて、akを演算する。その結果、
a（鉄相当）= 128349
a（アクリル樹脂相当）= 8885
が得られた。次いで、ステップ２０３において、物質鉄及びアクリル樹脂相当のＸ線吸収係数μikを数９の式を用いて演算する。次いで、ステップ１０７において、上述の長さLk、図５のＸ線連続スペクトルのエネルギーEi、ステップ２０３にて得られたＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、１回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図３６の（Ｂ）に示すＣＴ画像を再構成する。同様にして、２回目補正、３回目補正、４回目補正、５回目補正を行うと、図３６の（Ｃ）、（Ｄ）、（Ｅ）、（Ｆ）に示すＣＴ画像が得られる。計算回数は増加するが、画像アーチファクトが徐々に減少する。
【００７７】
第２のシミュレーション例は３つの物質つまり鉄、アルミニウム、アクリル樹脂を用いた図１８に示す複数異種ピンモデルに適用する。但し、物質鉄、アルミニウム、アクリル樹脂のＸ線吸収係数μikは未知とする。
【００７８】
第２のシミュレーション例においては、ステップ２０１において、図５のＸ線連続スペクトルのエネルギーEiを記憶する。そして、ステップ１０３において、これらに基づいて単一スペクトル投影データgを生成し、ステップ１０５において、単一スペクトル投影データgから図３７の（Ａ）に示すＣＴ画像を再構成する。
【００７９】
次に、ＣＴ画像に１回目補正を行う。つまり、ステップ１０６において、ＣＴ画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のＸ線経路毎に長さLkを演算する。次いで、ステップ２０２Ａにおいて、akを演算する。その結果、
a（鉄相当）= 128349
a（アルミニウム相当）= 8885
a（アクリル樹脂相当）= 8885
が得られた。次いで、ステップ２０３において、物質鉄及びアクリル樹脂相当のＸ線吸収係数μikを数９の式を用いて演算する。次いで、ステップ１０７において、上述の長さLk、図５のＸ線連続スペクトルのエネルギーEi、ステップ２０３にて得られたＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、１回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図３７の（Ｂ）に示すＣＴ画像を再構成する。同様にして、２回目補正、３回目補正、４回目補正、５回目補正を行うと、図３７の（Ｃ）、（Ｄ）、（Ｅ）に示すＣＴ画像が得られる。計算回数は増加するが、やはり、画像アーチファクトが徐々に減少する。
【００８０】
第３のシミュレーション例は図３８に示す生体組織モデルに適用する。但し、生体組織のＸ線吸収係数μikは未知とする。
【００８１】
第３のシミュレーション例においては、ステップ２０１において、図５のＸ線連続スペクトルのエネルギーEiを記憶する。そして、ステップ１０３において、これらに基づいて単一スペクトル投影データgを生成し、ステップ１０５において、単一スペクトル投影データgから図３９の（Ａ）に示すＣＴ画像を再構成する。
【００８２】
次に、図３９の（Ａ）に示すＣＴ画像に１回目補正を行う。つまり、ステップ１０６において、ＣＴ画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のＸ線経路毎に長さLkを演算する。次いで、ステップ２０２Ａにおいて、akを演算する。その結果、
a（鉄相当）= 128349
a（アクリル樹脂相当）= 8885
が得られた。次いで、ステップ２０３において、物質鉄及びアクリル樹脂相当のＸ線吸収係数μikを数９の式を用いて演算する。次いで、ステップ１０７において、上述の長さLk、図５のＸ線連続スペクトルのエネルギーEi、ステップ２０３にて得られたＸ線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ１０９において、１回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’からＣＴ画像を再構成する。同様にして、複数回補正を行うと、図３９の（Ｂ）に示すＣＴ画像が得られる。計算回数は増加するが、やはり、画像アーチファクトが徐々に減少する。
【符号の説明】
【００８３】
１：Ｘ線発生器
２：Ｘ線検出器
３：対象撮影物
４：制御回路
Iin：入射時Ｘ線エネルギー
Iout：透過後Ｘ線エネルギー
μik：Ｘ線吸収係数
g：単一スペクトル投影データ
f：連続スペクトル投影データ
ht：補正量
【図１】
2011203160.tif 000019
【図２】
2011203160.tif 000020
【図８】
2011203160.tif 000021
【図１０】
2011203160.tif 000022
【図１６】
2011203160.tif 000023
【図２０】
2011203160.tif 000024
【図２４】
2011203160.tif 000025
【図２５】
2011203160.tif 000026
【図２６】
2011203160.tif 000027
【図３１】
2011203160.tif 000028
【図３８】
2011203160.tif 000029
【図４０】
2011203160.tif 000030
【図３】
2011203160.tif 000031
【図４】
2011203160.tif 000032
【図５】
2011203160.tif 000033
【図６】
2011203160.tif 000034
【図７】
2011203160.tif 000035
【図９】
2011203160.tif 000036
【図１１】
2011203160.tif 000037
【図１２】
2011203160.tif 000038
【図１３】
2011203160.tif 000039
【図１４】
2011203160.tif 000040
【図１５】
2011203160.tif 000041
【図１７】
2011203160.tif 000042
【図１８】
2011203160.tif 000043
【図１９】
2011203160.tif 000044
【図２１】
2011203160.tif 000045
【図２２】
2011203160.tif 000046
【図２３】
2011203160.tif 000047
【図２７】
2011203160.tif 000048
【図２８】
2011203160.tif 000049
【図２９】
2011203160.tif 000050
【図３０】
2011203160.tif 000051
【図３２】
2011203160.tif 000052
【図３３】
2011203160.tif 000053
【図３４】
2011203160.tif 000054
【図３５】
2011203160.tif 000055
【図３６】
2011203160.tif 000056
【図３７】
2011203160.tif 000057
【図３９】
2011203160.tif 000058
【図４１】
2011203160.tif 000059
【図４２】
2011203160.tif 000060
【図４３】
2011203160.tif 000061