X線CT画像再構成方法及びX線CT画像再構成プログラム
【課題】試料のX線吸収差により、画像再構成に生じる放射状パターンのアーチファクトを低減する。
【解決手段】X線連続スペクトルを記憶し、単一スペクトルのX線投影データgを生成し、補正量を初期化し、X線投影データを補正量h、tで補正して撮影物に含む2つの物質のX線吸収係数のエネルギー依存性を考慮して撮影物に単一スペクトルのX線を照射した場合に相当する補正X線投影データg’を得、単一波長用画像再構成アルゴリズムを用いて逆投影演算してCT画像を再構成し、各物質を区別する1つのしきい値TH1、TH2、…を用いてCT画像から各物質の領域を抽出しこの領域のX線経路の長さを演算し、発生X線連続スペクトルのエネルギー、各物質の領域の長さ及び各物質のX線吸収係数を用いて連続スペクトルのX線投影データfを演算し、単一ペクトルのX線投影データgと連続スペクトルのX線投影データfとの差で補正量h、tを演算する。
【解決手段】X線連続スペクトルを記憶し、単一スペクトルのX線投影データgを生成し、補正量を初期化し、X線投影データを補正量h、tで補正して撮影物に含む2つの物質のX線吸収係数のエネルギー依存性を考慮して撮影物に単一スペクトルのX線を照射した場合に相当する補正X線投影データg’を得、単一波長用画像再構成アルゴリズムを用いて逆投影演算してCT画像を再構成し、各物質を区別する1つのしきい値TH1、TH2、…を用いてCT画像から各物質の領域を抽出しこの領域のX線経路の長さを演算し、発生X線連続スペクトルのエネルギー、各物質の領域の長さ及び各物質のX線吸収係数を用いて連続スペクトルのX線投影データfを演算し、単一ペクトルのX線投影データgと連続スペクトルのX線投影データfとの差で補正量h、tを演算する。
【発明の詳細な説明】
【技術分野】
【0001】
本発明はX線CT(Computed Tomography)画像再構成方法及びX線CT画像再構成プログラムに関する。
【背景技術】
【0002】
一般に、X線CT画像再構成方法においては、撮影対象物の周囲を多数の方向から通常の波長分布を有するX線を投影し、X線の透過エネルギーとしての投影データを取得し、この投影データをフィルタ補正逆投影法、逐次近似法等の単一波長用画像再構成アルゴリズムを用いて逆投影して画像(CT画像)を再構成する。
【0003】
上述のX線CT画像再構成方法において、撮影対象物にX線吸収係数が極端に大きい物質領域たとえば金属が存在すると、その物質領域を波長分布を有する照射X線が透過することによりその線質が変化するビームハードニング現象により、画像アーチファクトと呼ばれる放射状ノイズパターンがCT画像に発生する。この画像アーチファクトは医療現場、非破壊検査現場等で再構成画像の判別の妨げとなる。
【0004】
次に、ビームハードニング現象について、図40、図41、図42、図43を参照して説明する。
【0005】
図40は一般的なX線CT画像再構成装置を示す図である。
【0006】
図40において、1はX線発生器、2はX線検出器、3は撮影対象物、4は、X線発生器1、X線検出器2、撮影対象物3を制御する制御回路である。制御回路4はコンピュータで構成され、CPU、ROM、RAM等を含む。
【0007】
図41は図40の撮影対象物3の物質のX線吸収係数μを示す。
【0008】
図41に示すように、X線吸収係数μの分布が物質毎に異なるだけでなく、同一物質でもX線波長によってX線吸収係数μが異なる、つまり、エネルギー依存性を有する。
【0009】
図42は図40のCT画像再構成装置のX線エネルギースペクトルを示し、(A)はX線発生器1の発生X線の連続エネルギースペクトル(以下、単に、連続スペクトル)を示し、(B)は大きさL、X線吸収係数μ(E)の撮影対象物3を透過後のX線検出器2の検出X線の連続スペクトルを示す。
【0010】
図42に示すように、撮影対象物3の透過前後でX線の連続スペクトルは異なる。
【0011】
ところで、図40のX線CT画像再構成装置におけるCT画像再構成演算においては、X線吸収量のエネルギー依存性を考慮していない。つまり、μ一定の単一エネルギースペクトル(以下単に、単一スペクトル)のX線エネルギー照射つまり単一波長照射を前提としている。この理由は、X線連続スペクトルは装置、照射するX線等の種類、撮影条件等に依存するので、これに対応することが困難であるためである。しかし、実際のCT撮影は通常のX線のために波長分布を有する。この結果、単一波長用CT画像再構成時の演算によって矛盾が発生し、図43に示すような画像アーチファクトが発生する。
【0012】
画像アーチファクトを低減する第1の従来のX線CT画像再生方法は投影データのシノグラムから高X線吸収係数の物質領域たとえば金属領域を抽出して仮想投影によりその物質領域の投影データを改変して修正投影データを演算し、この修正投影データを逆投影して再構成画像(CT画像)を再構成している(参照:特許文献1)。
【0013】
上述の第1の従来のX線CT画像構成方法においては、厳しい画像アーチファクトを低減することができるが、軟組織部分の抽出が難しく投影データの修正が不充分であった。
【0014】
画像アーチファクトを低減する第2の従来のX線CT画像構成方法は、投影条件を変更して複数回CT撮影を行って複数の投影データを得てX線連続スペクトルの種類を限定し、X線吸収量を求めている(参照:特許文献2、非特許文献1)。
【先行技術文献】
【特許文献】
【0015】
【特許文献1】特開2006−167161号公報
【特許文献2】特開2008−241376号公報
【非特許文献】
【0016】
【非特許文献1】李根旭、小関道彦、木村仁、伊能教夫:X線CT画像の画質向上に関する研究(X線スペクトルを考慮した画像再構成手法の提案)、計測自動制御学会 第24回センシングフォーラム資料、pp.107-111,2007
【発明の概要】
【発明が解決しようとする課題】
【0017】
しかしながら、上述の第2の従来のX線CT画像再構成方法においては、複数のX線連続スペクトルを必要とするので、CT撮影回数が増加すると共に、撮影条件を変更することによる異なるX線連続スペクトルを得ることが難しいという課題がある。
【課題を解決するための手段】
【0018】
上述の課題を解決するために本発明は、X線CT画像再構成方法において、発生X線連続エネルギースペクトルを既知情報として記憶し、撮影対象物の単一エネルギースペクトルのX線投影データを生成し、補正量を初期化し、単一エネルギースペクトルのX線投影データを補正量で補正して撮影対象物に含まれる少なくとも2つの物質のX線吸収係数のエネルギー依存性を考慮して撮影対象物に単一エネルギースペクトルのX線を照射した場合に相当する補正X線投影データを得、補正X線投影データを単一波長用画像再構成アルゴリズムを用いて逆投影演算してCT画像を再構成し、各物質を区別するための少なくとも1つのしきい値を用いてCT画像から各物質の領域を抽出してこれらの領域のX線経路の長さを演算し、発生X線連続エネルギースペクトルのエネルギー、各物質の領域の長さ及び各物質のX線吸収係数を用いて連続エネルギースペクトルのX線投影データを演算し、単一エネルギースペクトルのX線投影データと連続エネルギースペクトルのX線投影データとの差により補正量を演算するものである。すなわち、従来の単一波長用の画像再構成アルゴリズムを用いて得られるCT画像(元CT画像)に対して、推定される撮影対象物に対して単一エネルギースペクトルのX線を照射したときX線投影データを推定する。この推定には、元CT画像に少なくとも1つのしきい値を用いて少なくとも2つの物質の領域を抽出し、各X線経路の長さを演算する。次いで、既知情報のX線連続エネルギースペクトルのエネルギー、物質の領域の長さ及びX線吸収係数を用いて実際に照射された連続エネルギースペクトルのX線投影データを補正演算し、単一エネルギースペクトルのX線投影データに近い値を得るようにしたものである。このようにして、少なくとも1つのしきい値で区分された少なくとも2つの物質が及ぼすビームハードニングの度合を個々に換算して補正量で表わす。
【0019】
また、上述の補正から補正量演算を繰返して補正量が所定範囲内に収束させる。この場合、しきい値が未知であれば、CT画像のヒストグラムによりしきい値を推定してしきい値の数を増加させる。このようにして、補正量が更新され、画像アーチファクトが徐々に減少することで物質の領域のX線経路の長さが正確となる。
【発明の効果】
【0020】
本発明によれば、単一エネルギースペクトルのX線投影データは変更されないので、CT撮影は1回で済み、しかも、補正を繰返すことによって画像アーチファクトを低減することができる。
【図面の簡単な説明】
【0021】
【図1】本発明に係るX線CT画像再構成方法の第1の実施の形態を説明するためのフローチャートである。
【図2】図1のステップ101のX線連続スペクトルの例を示すグラフである。
【図3】図1のステップ106のしきい値を説明する図であって、(A)はCT画像、(B)は画素値ヒストグラムである。
【図4】図1のステップ108の補正量を説明するためのグラフである。
【図5】X線発生器の発生X線の連続スペクトルを示す表である。
【図6】生体骨のX線吸収係数を示す表である。
【図7】図1の第1のシミュレーション例に用いられる撮影対象物としての骨円柱モデルを示す図である。
【図8】図7の補正前の骨円柱モデルを示し、(A)はCT画像、(B)は(A)センタライン上のCT値を示すグラフである。
【図9】図8の(A)のCT画像から物質の領域抽出を示すCT画像である。
【図10】図7の1回目補正の骨円柱モデルを示し、(A)はCT画像、(B)は(A)のセンタライン上のCT値を示すグラフである。
【図11】鉄のX線吸収係数を示す表である。
【図12】アクリル樹脂のX線吸収係数を示す表である。
【図13】図1の第2のシミュレーション例に用いられる撮影対象物としての複数同種ピンモデルを示す図であって、(A)は斜視像、(B)は上面像である。
【図14】図13の補正前の複数同種ピンモデルを示し、(A)はCT画像、(B)は(A)センタライン上のCT値を示すグラフである。
【図15】図14の(A)のCT画像から物質の領域抽出を示し、(A)は鉄の領域抽出を示すCT画像、(B)はアクリル樹脂の領域抽出を示すCT画像である。
【図16】図13の1回目補正の複数同種ピンモデルを示し、(A)はCT画像、(B)は(A)のセンタライン上のCT値を示すグラフである。
【図17】アルミニウムのX線吸収係数を示す表である。
【図18】図1の第3のシミュレーション例に用いられる撮影対象物としての複数異種ピンモデルを示す上面像である。
【図19】図18の補正前の複数異種ピンモデルを示し、(A)はCT画像、(B)は(A)センタライン上のCT値を示すグラフである。
【図20】図19の画素値ヒストグラムを示すグラフである。
【図21】図19の(A)のCT画像から物質の領域抽出を示し、(A)は鉄の領域抽出を示すCT画像、(B)はアルミニウムの領域抽出を示すCT画像、(C)はアクリル樹脂の領域抽出を示すCT画像である。
【図22】図18の1回目補正の複数異種ピンモデルを示し、(A)はCT画像、(B)は(A)の画素値ヒストグラムを示すグラフである。
【図23】図18の2回目補正の複数異種ピンモデルを示し、(A)はCT画像、(B)は(A)の画素値ヒストグラムを示すグラフである。
【図24】図18の3回目補正の複数異種ピンモデルを示し、(A)はCT画像、(B)は(A)の画素値ヒストグラムを示すグラフである。
【図25】図1の第3のシミュレーションの補正量の収束を示すグラフである。
【図26】図1の第1の実験例を示し、(A)は補正前のCT画像、(B)は補正後のCT画像である。
【図27】図1の第2の実験例を示し、(A)は補正前のCT画像、(B)は補正後のCT画像である。
【図28】図1の第3の実験例を示し、(A)は鉄釘―米粒モデルの斜視像、(B)は撮影時のX線スペクトルを示すグラフである。
【図29】図1の第3の実験例を示し、(A)は補正前のCT画像、(B)は補正後のCT画像である。
【図30】第3の実験例において鉄のみ抽出による補正後のCT画像である。
【図31】本発明に係るX線CT画像再構成方法の第2の実施の形態を説明するためのフローチャートである。
【図32】図31のステップ202、203における物質鉄、アルミニウム、アクリル樹脂のX線吸収係数の例を示すグラフである。
【図33】図31の第1のシミュレーション例の1回目補正後のCT画像である。
【図34】図31の第2のシミュレーション例を示し、(A)は1回目補正後のCT画像、(B)は補正後のCT画像である。
【図35】本発明に係るX線CT画像再構成方法の第3の実施の形態を説明するためのフローチャートである。
【図36】図35の第1のシミュレーション例を示し、(A)、(B)、(C)、(D)、(E)、(F)は、補正前、1回目補正後、2回目補正後、3回目補正後、4回目補正後、5回目補正後のCT画像である。
【図37】図35の第2のシミュレーション例を示し、(A)、(B)、(C)、(D)、(E)は、補正前、1回目補正後、2回目補正後、3回目補正後、4回目補正後のCT画像である。
【図38】図35の第1のシミュレーション例に用いられる撮影対象物としての生体組織モデルを示す図である。
【図39】図38のCT画像を示し、(A)は補正前のCT画像、(B)は補正後のCT画像である。
【図40】一般的なX線CT画像再構成装置を示す図である。
【図41】図40の撮影対象物の物質のX線吸収係数を示すグラフである。
【図42】図40のX線連続スペクトルを示すグラフである。
【図43】画像アーチファクトの例を示すCT画像である。
【発明を実施するための形態】
【0022】
図1は本発明に係るX線CT画像再構成方法の第1の実施の形態を説明するためのフローチャートである。図1のフローチャートは図40のX線CT画像再構成装置によって実行される。
【0023】
始めに、ステップ101では、X線発生器1のX線連続スペクトルを予め記憶する。このX線発生器1のX線管の種類、撮影条件によって定まるX線連続スペクトルは専用の機器で予め計測できる。尚、X線連続スペクトルは実際には離散値としてEi及びその光子数n(Ei)で記憶される。
【0024】
次に、ステップ102では、後述の補正量ht及び補正回数tを初期化する。つまり、
h0←0
t←0
ここで、h0は投影データの数だけ存在するので、h0はその総称を表している。
【0025】
次にステップ103では、X線検出器2の検出信号により単一スペクトルつまり波長iでの単一スペクトル投影データ(X線吸収量)gを生成する。尚、波長iはX線管電圧から計算される代表値を採用する。
【0026】
ステップ103におけるX線吸収量gを理論的に考察すると、X線吸収量gは数1の式で表される。
【数1】
但し、Iin:物質kの入射時のエネルギー
Iout:物質kの透過後のエネルギー
Lk:物質kのX線経路の長さ
μik:波長i及び物質kに依存するX線吸収係数
X線が透過する物質がn種類の物質で構成されていれば、数1の式は数2の式となる
【数2】
数2の式を対数で表すと、数3の式となる。
【数3】
【0027】
次に、ステップ104では、単一スペクトル投影データgに補正量htを減算して補正投影データg’を求める。
g’= g - ht
尚、ここでも、g’は投影データの数だけ存在するので、その総称を表している。
【0028】
次に、ステップ105では、補正投影データg’をフィルタ補正逆投影法等により逆投影してCT画像を再構成する。この場合の逆投影法は、従来の単一スペクトルのX線投影データを再構成するための単一波長用画像再構成アルゴリズムを用いる。
【0029】
次に、ステップ106では、既知のしきい値TH1、TH2、…を用いてCT画像から物質kのX線経路の領域の長さLkを演算する。たとえば、図3の(A)に示すCT画像が図3の(B)に示すごとく3種類の物質、つまり、鉄、アルミニウム、アクリル樹脂を含んでいれば、ステップ106において、TH1を画素値14000、TH2を画素値600に予め設定しておく。尚、対象撮影物3が複数の未知の物質を含んでいる場合には、ステップ106において、図3の(B)の画素値ヒストグラムを演算し、画素値の集中箇所、たとえば、図3の(B)における画素値-900、1800、8000、20000の集中箇所を判別し、その間にしきい値を推定して設定する。尚、このしきい値推定アルゴリズムは容易に得ることができる。
【0030】
次に、ステップ107では、各単一スペクトル投影データgに対応する連続スペクトル投影データfを演算する。
【0031】
連続スペクトル投影データfはステップ101にて記憶されたX線発生器1のX線連続スペクトルのエネルギーEi、物質kのX線吸収係数μik及びステップ106にて演算された物質kのX線経路の長さLkに基づいて演算される。つまり、物質kの入射時のX線エネルギーIinは数4の式で表わされる。
【数4】
また、物質kの透過後のX線エネルギーIoutは数5の式で表わされる。
【数5】
数4、数5の式から数6の式が得られる。
【数6】
数6の式を対数で表わすと、連続スペクトル投影データ(X線吸収量)fは数7の式で表わされる。
【数7】
【0032】
次に、ステップ108では、
t ← t + 1
とする。
【0033】
次に、ステップ109では、t回目の補正量htを演算する。すなわち、
ht = g - f
【0034】
すなわち、図4に示すように、物質kのX線経路の長さLkが同一であっても、単一スペクトル投影データ(X線吸収量)gと連続スペクトル投影データ(X線吸収量)fとでは値が補正量htだけ異なっている。この補正量htだけ異なっている分、この補正量htが画像アーチファクトとして現われている。ステップ106における画像アーチファクトが存在する投影データg’からの物質kのX線経路の長さLkにも僅差が存在する。従って、ステップ104において単一スペクトル投影データgをこの補正量htで減算して補正することにより連続スペクトル投影データfに近づけることができ、長さLkも安定的収束し、これに伴い、画像アーチファクト低減効果が生ずる。
【0035】
ステップ110はステップ104〜109のフローを補正量htが所定範囲に収束して|ht-ht-1|がε未満になるまで繰返す。|ht-ht-1|がε未満となったときにステップ111にてこのルーチンは終了する。尚、ここでも、|ht-ht-1|は投影データの数だけ存在し、その総和を意味する。
【0036】
図1のフローチャートで重要なのは、1回のCT撮影で得られた単一スペクトル投影データgは変更されず、その補正量htのみを更新し、この更新によって画像アーチファクトを徐々に低減することで物質kのX線経路の長さLkの演算も正確になっていくことである。
【0037】
図1のフローチャートによるシミュレーション例を説明する。つまり、この場合、制御回路(コンピュータ)4上でCT撮影する仮想投影手法を用い、X線発生器1のX線連続スペクトルは図5に示す9個の離散値で近似する。
【0038】
第1のシミュレーション例は図6に示すX線吸収係数μikを有する物質つまり生体骨を用いた図7に示す骨円柱モデルに適用する。つまり、物質k及びそのX線吸収係数μikは既知とする。この場合、骨円柱モデルは均一であり、かつ形状が単純であるので、ビームハードニングの影響が小さく、従って、画像アーチファクトが少ないと考えられるケースである。
【0039】
第1のシミュレーション例においては、ステップ101において、図5のX線連続スペクトル及び図6のX線吸収係数μikを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。この場合、図8の(A)に示すCT画像が得られ、一見すると画像アーチファクトは出現していないが、図8の(B)に示すように、CT画像のセンタライン上のCT値をプロットすると、円柱中央部の画素で輝度が上昇し、円柱周辺部の画素で輝度が低下していることが分かる。つまり、円柱周辺部の画素のCT値が真値1240より大きくなっている。このように、単純な形状でもビームハードニングの影響で適切なX線吸収量が得られず、CT値の定量性が喪失してしまう。
【0040】
次に、図8の(A)のCT画像に1回目補正を行う。つまり、ステップ106において、図8の(A)のCT画像の画素に対して既知のしきい値TH1たとえば500を用いて図9に示す物質生体骨の領域を抽出する。図9に示すごとく、物質生体骨の領域は正確に抽出されている。この抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ107において、その長さLk、図5のX線連続スペクトルのエネルギーEi、図6のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データから図10の(A)に示すCT画像を再構成する。
【0041】
図10の(A)の1回目補正のCT画像についてセンタライン上の画素をプロットすると、図10の(B)に示すごとく、円柱周辺部の画素の輝度が上昇して円柱内部で均一となり、50keVのCT真値1240に近づいていることが分かる。
【0042】
第2のシミュレーション例は図11、図12に示すX線吸収係数μikを有する物質つまり鉄、アクリル樹脂を用いた図13に示す複数同種ピンモデルに適用する。類似モデルとしては、電子部品、半導体の非破壊検査の際に、金属部品が樹脂中に埋め込められていると、金属部品の周囲に画像アーチファクトが発生して金属部品の周囲の樹脂部分の形状、組成が判別できなくなる場合がある。また、人体にインプラント金属が埋め込められていると、インプラント金属の周囲に画像アーチファクトが発生してインプラント金属の周囲の筋肉部分の形状、組成が判別できなくなる場合がある。尚、この場合も、物質k及びそのX線吸収係数μikは既知とする。この複数同種ピンモデルは形状が少し複雑であるので、ビームハードニングの影響が大きく、従って、画像アーチファクトがある程度大きいと考えられるケースである。
【0043】
第2のシミュレーション例においては、ステップ101において、図5のX線連続スペクトルのエネルギーEi及び図11、図12のX線吸収係数μikを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。この場合、図14の(A)に示すCT画像が得られ、強い画像アーチファクトが出現している。図14の(B)に示すように、CT画像のセンタライン上のCT値をプロットすると、ピン中央部の画素で輝度が急上昇し、ピン周辺部の画素で輝度が急低下していることが分かる。つまり、ピン周辺部の画素のCT値が真値20000より大きくなっている。このように、複雑な形状では強いビームハードニングの影響で適切なX線吸収量が得られず、CT値の定量性が喪失してしまう。
【0044】
次に、図14の(A)のCT画像に1回目補正を行う。つまり、ステップ106において、図14の(A)のCT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて図15に示す物質鉄及びアクリル樹脂の各領域を抽出する。この場合、図15の(A)に示すごとく、物質鉄の領域は正確に抽出されているが、図15の(B)に示すごとく、物質アクリル樹脂の領域は物質鉄に起因するアーチファクトの影響で不正確である。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ107において、これらの長さLk、図5のX線連続スペクトルのエネルギーEi、図11、図12のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図16の(A)に示すCT画像を再構成する。
【0045】
図16の(A)の1回目補正のCT画像についてセンタライン上画素をプロットすると、図16の(B)に示すごとく、ピン周辺部の画素の輝度が上昇してピン内部で均一となり、50keVのCT真値20000に近づいていることが分かる。また、同時に、物質アクリル樹脂の領域の画素の輝度も均一となり、50keVのCT真値-680に近づいていることが分かる。
【0046】
第3のシミュレーション例は図11、図12、図17に示すX線吸収係数μikを有する物質つまり鉄、アクリル樹脂、アルミニウムを用いた図18に示す複数異種ピンモデルに適用する。尚、この場合も、物質k及びそのX線吸収係数μikは既知とする。この複数異種ピンモデルは複数の高X線吸収係数を有する物質を有するので、ビームハードニングの影響が非常に大きく、従って、画像アーチファクトが大きいと考えられるケースである。
【0047】
第3のシミュレーション例においては、ステップ101において、図5のX線連続スペクトルのエネルギーEi及び図11、図12、図17のX線吸収係数μikを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。この場合、図19の(A)に示すCT画像が得られ、物質鉄に起因する強い画像アーチファクトが出現している。図19の(B)に示すように、CT画像のセンタライン上のCT値をプロットすると、ピン中央部の画素で輝度が急上昇し、ピン周辺部の画素で輝度が急低下していることが分かる。つまり、ピン周辺部の画素のCT値が真値20000あるいは1800より大きくなっている。このときの画素ヒストグラムを図20に示すと、物質鉄の分布は真値20000から離れた23000程度の箇所にある。従って、図1のステップ106におけるしきい値TH1を8000と推定して設定し、物質鉄の領域を抽出できる。尚、この場合、物質鉄が未知であっても、1つのしきい値TH1が推定され設定される。この結果、図21の(A)に示すごとく、物質鉄の領域を正確に抽出できるが、図21の(B)、(C)に示すごとく、物質アルミニウム、アクリル樹脂の各領域は物質鉄に起因するアーチファクトの影響で正確には抽出できない。
【0048】
次に、図19の(A)のCT画像に1回目補正を行う。つまり、上述の抽出領域のX線経路毎に長さLkを演算する。次いで、これらの長さLk、図5のX線連続スペクトルのエネルギーEi、図11、図12、図17のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図22の(A)に示すCT画像を再構成する。図22の(A)のCT画像においては、画像アーチファクトが残り、補正は不充分である。この理由は、物質鉄に起因する画像アーチファクトの画素値が物質アルミニウムの画素値とほぼ同一であるために、画像アーチファクトの形状が物質アルミニウムが存在するものとして再構成されたものと考えられる。このときの画素ヒストグラムを図22の(B)に示すと、物質アルミニウム及びアクリル樹脂のピーク値がこれらの真値-680及び1800に現われている。従って、図1のステップ106におけるしきい値TH2を-680と1800との間たとえば1200と推定して設定すれば、物質アルミニウム、アクリル樹脂の各領域は正確に抽出できる。尚、この場合も、物質アルミニウム、アクリル樹脂が未知であっても、しきい値TH2が推定され設定される。従って、撮影対象物に含まれるn個の物質が未知である場合には、図1のステップ104〜110を繰返すことによりn−1個のしきい値が推定され設定されることになる。この結果、未知の物質の領域も正確に抽出できる。
【0049】
次に、図22の(A)のCT画像に2回目補正を行う。つまり、ステップ106において、図22の(A)のCT画像の画素に対して推定のしきい値TH1たとえば8000及びしきい値TH2たとえば1200を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。次いで、ステップ107において、これらの長さLk、図5のX線連続スペクトルのエネルギーEi、図11、図12、図17のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、2回目の補正量h2を演算し、補正量h2で単一スペクトル投影データgを補正し、補正された投影データg’から図23の(A)に示すCT画像を再構成する。
【0050】
次に、図23の(A)のCT画像に3回目補正を行う。つまり、ステップ106において、図23の(A)のCT画像の画素に対して推定のしきい値TH1たとえば8000及びしきい値TH2たとえば1200を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。次いで、ステップ107において、これらの長さLk、図5のX線連続スペクトルのエネルギーEi、図11、図12、図17のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、3回目の補正量h3を演算し、補正量h3で単一スペクトル投影データgを補正し、補正された投影データg’から図23の(A)に示すCT画像を再構成する。
【0051】
図24の(A)の3回目補正のCT画像についてセンタライン上画素をプロットすると、図24の(B)に示すごとく、ピン周辺部の画素の輝度が上昇してピン内部で均一となり、50keVのCT真値20000もしくは1800に近づいていることが分かる。また、同時に、物質アクリル樹脂の領域の画素の輝度も均一となり、50keVのCT真値-680に近づいていることが分かる。
【0052】
第3のシミュレーション例においては、CT画像をたとえ4回目補正、5回目補正、…を行ってもCT画像はほとんど変わらない。つまり、補正回数が増加するにつれて、画像アーチファクトが低減し、物質鉄、アルミニウム、アクリル樹脂の領域がほぼ正確に抽出できる。従って、図25に示すごとく、補正回数が3回を超えて総和|ht-ht-1|がε以下となると、各物質がその真値に近づいており、CT画像の誤差は低減し、補正の効果が向上していることが確認できる。このとき、補正量と前回補正量との差|ht-ht-1|も0に近づき、補正量が一定値に収束していることが分かる。
【0053】
次に、図1のフローチャートによる実験例を説明する。
【0054】
図26は上述の第2のシミュレーション例に対応して行った第1の実験例である。撮影時の発生X線の連続スペクトルは不明であったので、図5の連続スペクトルの離散値で近似させた。また、物質鉄、アクリル樹脂のX線吸収係数μikは図11、図12に示す値を既知として用いた。さらに、物質鉄、アクリル樹脂の各領域の抽出のためのしきい値に適当な値を既知として用いた。この結果、補正前のCT画像は図26の(A)に示すごとくであったが、図1のフローチャートに基づく補正を行った結果、図26の(B)に示すCT画像が得られた。つまり、画像アーチファクトは大幅に減少し、補正の効果が認められた。
【0055】
図27は鉄の十字架形状モデルに対して行った第2の実験例である。この場合も、撮影時の発生X線の連続スペクトルは不明であったので、図5の連続スペクトルの離散値で近似させた。また、物質鉄、アクリル樹脂のX線吸収係数μikは図11、図12に示す値を既知として用いた。さらに、物質鉄、アクリル樹脂の各領域の抽出のためのしきい値に適当な値を既知として用いた。この結果、補正前のCT画像は図27の(A)に示すごとくであったが、図1のフローチャートに基づく補正を行った結果、図27の(B)に示すCT画像が得られた。やはり、画像アーチファクトは大幅に減少し、補正の効果が認められた。
【0056】
図28の(A)に示す9本の鉄釘を米粒の山へ埋込んだモデルに対して行ったものである。撮影時の管電圧130kVの発生X線の連続スペクトルは、図28の(B)の連続スペクトルの離散値で近似させた。また、物質鉄、米粒のX線吸収係数μikは図11、図12に示す値を既知として用いた。つまり、米粒はアクリル樹脂とみなした。さらに、物質鉄、米粒の各領域の抽出のためのしきい値に13500及び750を既知として用いた。この結果、補正前のCT画像は図29の(A)に示すごとくであった。つまり、鉄釘と鉄釘とを結ぶラインで強い画像アーチファクトが発生している。これに対し、図1のフローチャートに基づく補正を行った結果、図29の(B)に示すCT画像が得られた。やはり、画像アーチファクトは大幅に減少し、補正の効果が認められた。
【0057】
尚、第3の実験例において、しきい値13500だけ設定して鉄釘(鉄)のみの領域を抽出した場合には、図30に示すごとく、米粒が正確に抽出できなかった。この理由は、連続スペクトルにより、鉄釘(鉄)だけでなく、米粒部分でもX線吸収量の過吸収が発生しているものと考えられるからである。
【0058】
このように、実験例によれば、本発明はCT画像を用いた個体別モデリング、リバースエンジニアリング非破壊検査で有効であることが分かる。
【0059】
図31は本発明に係るX線CT画像再構成方法の第2の実施の形態を説明するためのフローチャートであって、対象撮影物3を構成する物質kの数が少なくかつこれら物質kのX線吸収係数μikが未知の場合に適用される。
【0060】
図31においては、図1において、ステップ101をステップ201に置換し、ステップ202、203を追加したものである。すなわち、ステップ201では、図5のX線連続スペクトルのエネルギーEiのみを記憶し、ステップ202、203において、X線連続スペクトルのエネルギーEi及び物質kのX線経路の長さLkに基づいて物質kのX線吸収係数μikを演算する。
【0061】
ステップ202、203の根拠を説明する。X線検出器2のある位置の投影データDは数8の式で表わされる。
【数8】
ここで、X線連続スペクトルのエネルギーEiが既知であるので、数8の式の投影データはX線吸収係数μik及びステップ106で定まる長さLkの関数となる。他方、X線吸収係数μikは数9の式で表わされる。
【数9】
数9の式の妥当性は、図32に示すように、物質鉄、アルミニウム、アクリル樹脂に対してak=134000、7000、200としたときに、数9の式により求められたX線吸収係数μikは実際のX線吸収係数μikにほぼ一致することから確認できる。
数9の式で数8の式を表わすと、数10の式となる。
【数10】
数10の式において、物質kの種類が未知数n個のn次方程式としてakを解くことができる。実際に検出されるデータはX線検出器2の列数たとえば512×撮影の回転方向数たとえば800個分入手できるので、物質kの種類数nはどんなに多くとも解くことができる。つまり、複数の方向から照射したX線ビーム毎に投影データD’が得られるので、数10の式がその数だけ成立する。この場合、数10の式の数はn以上であればakを求めることができる。但し、領域の長さLkは画像アーチファクトによる誤差を有するので、akは最小2乗法で求める。尚、このとき、求められたakに基づき、物質kを推定できる。
【0062】
従って、ステップ202では、数10の式のakを最小2乗法により演算する。次いで、ステップ203において、物質kのX線吸収係数μikを数9の式により演算する。
【0063】
図31のフローチャートによるシミュレーション例を説明する。つまり、この場合、制御回路(コンピュータ)4上でCT撮影する仮想投影手法を用い、X線発生器1のX線連続スペクトルは図5に示す9個の離散値で近似する。
【0064】
第1のシミュレーション例は2つの物質つまり鉄、アクリル樹脂を用いた図13に示す複数同種ピンモデルに適用する。但し、物質鉄、アクリル樹脂のX線吸収係数μikは未知とする。
【0065】
第1のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。
【0066】
次に、CT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202において、数10の式のakを最小2乗法により演算する。その結果、
a1 = 124651
a2 = 41
が得られた。尚、a(鉄)、a(アクリル樹脂)の真値は130000、200であるので、上述の結果値は真値に近いことが分かる。従って、a1 = 124651は鉄に相当、a2= 41はアクリル樹脂に相当すると推定する。次いで、ステップ203において、物質鉄及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図33に示すCT画像を再構成する。このように、X線吸収係数μikが未知である場合にも、画像アーチファクトが大幅に減少する。
【0067】
第2のシミュレーション例は3つの物質つまり鉄、アルミニウム及びアクリル樹脂を用いた図18に示す複数異種ピンモデルに適用する。但し、物質鉄、アルミニウム及びアクリル樹脂のX線吸収係数μikは未知とする。
【0068】
第2のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。
【0069】
次に、CT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202において、数10の式のakを最小2乗法により演算する。その結果、
a1 = 127682
a2 = -1702
a3 = 416
が得られた。この場合、画像アーチファクトの影響によりアルミニウム相当の領域抽出が不適切であるために、a(アルミニウム)の真値7000から離れてしまった。また、X線吸収係数は負とならないので、
a1 = 127682
a2 = 0
a3 = 416
とする。次いで、ステップ203において、物質鉄、アルミニウム及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図34の(A)に示すCT画像を再構成する。このように、1回目補正により、不充分であるが、画像アーチファクトが減少する。
【0070】
次に、図34の(A)のCT画像に2回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202において、数10の式のakを最小2乗法により演算する。その結果、
a1 = 123476
a2 = 6962
a3 = 21
が得られた。尚、a(鉄)、a(アルミニウム)、a(アクリル樹脂)の真値は130000、7000、200であるので、上述の結果値は真値に近いことが分かる。従って、a1 = 123476は鉄、a2= 6962はアルミニウム、a3 = 21はアクリル樹脂に相当すると推定する。次いで、ステップ203において、物質鉄、アルミニウム及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図34の(B)に示すCT画像を再構成する。このように、2回目補正により、画像アーチファクトが大幅に減少する。
【0071】
図35は本発明に係るX線CT画像再構成方法の第3の実施の形態を説明するためのフローチャートであって、対象撮影物3を構成する物質kの数が多くかつこれら物質kのX線吸収係数μikが未知の場合に適用される。
【0072】
図35においては、図31において、ステップ202をステップ202Aに置換したものである。すなわち、ステップ201Aでは、図31のフローチャートにおけるakを
ak = 9.873・THn + 9873
= 9.873・(THn+ 1000)
但し、THn(n=1,2,…)は暫定的なしきい値、
により演算する。ここで、値9.873、9873は、空気のCT値が-1000である場合にa(空気)=0、鉄のCT値12572の場合にa(鉄)=134000となるように線形に規定したものである。ここで、鉄のCT値12572は実際のCT画像の最大CT値である。但し、9.873、9873は与えられた2点が線形となるように、他の値αk、βkとなし得る。この場合も、求められたakに基づき、物質kを推定できる。
【0073】
図35のフローチャートによるシミュレーション例を説明する。つまり、この場合、制御回路(コンピュータ)4上でCT撮影する仮想投影手法を用い、X線発生器1のX線連続スペクトルは図5に示す9個の離散値で近似する。
【0074】
第1のシミュレーション例は2つの物質つまり鉄、アクリル樹脂を用いた図13に示す複数同種ピンモデルに適用する。但し、物質鉄、アクリル樹脂のX線吸収係数μikは未知とする。
【0075】
第1のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgから図36の(A)に示すCT画像を再構成する。
【0076】
次に、CT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202Aにおいて、akを演算する。その結果、
a(鉄相当)= 128349
a(アクリル樹脂相当)= 8885
が得られた。次いで、ステップ203において、物質鉄及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図36の(B)に示すCT画像を再構成する。同様にして、2回目補正、3回目補正、4回目補正、5回目補正を行うと、図36の(C)、(D)、(E)、(F)に示すCT画像が得られる。計算回数は増加するが、画像アーチファクトが徐々に減少する。
【0077】
第2のシミュレーション例は3つの物質つまり鉄、アルミニウム、アクリル樹脂を用いた図18に示す複数異種ピンモデルに適用する。但し、物質鉄、アルミニウム、アクリル樹脂のX線吸収係数μikは未知とする。
【0078】
第2のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgから図37の(A)に示すCT画像を再構成する。
【0079】
次に、CT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202Aにおいて、akを演算する。その結果、
a(鉄相当)= 128349
a(アルミニウム相当)= 8885
a(アクリル樹脂相当)= 8885
が得られた。次いで、ステップ203において、物質鉄及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図37の(B)に示すCT画像を再構成する。同様にして、2回目補正、3回目補正、4回目補正、5回目補正を行うと、図37の(C)、(D)、(E)に示すCT画像が得られる。計算回数は増加するが、やはり、画像アーチファクトが徐々に減少する。
【0080】
第3のシミュレーション例は図38に示す生体組織モデルに適用する。但し、生体組織のX線吸収係数μikは未知とする。
【0081】
第3のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgから図39の(A)に示すCT画像を再構成する。
【0082】
次に、図39の(A)に示すCT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202Aにおいて、akを演算する。その結果、
a(鉄相当)= 128349
a(アクリル樹脂相当)= 8885
が得られた。次いで、ステップ203において、物質鉄及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’からCT画像を再構成する。同様にして、複数回補正を行うと、図39の(B)に示すCT画像が得られる。計算回数は増加するが、やはり、画像アーチファクトが徐々に減少する。
【符号の説明】
【0083】
1:X線発生器
2:X線検出器
3:対象撮影物
4:制御回路
Iin:入射時X線エネルギー
Iout:透過後X線エネルギー
μik:X線吸収係数
g:単一スペクトル投影データ
f:連続スペクトル投影データ
ht:補正量
【技術分野】
【0001】
本発明はX線CT(Computed Tomography)画像再構成方法及びX線CT画像再構成プログラムに関する。
【背景技術】
【0002】
一般に、X線CT画像再構成方法においては、撮影対象物の周囲を多数の方向から通常の波長分布を有するX線を投影し、X線の透過エネルギーとしての投影データを取得し、この投影データをフィルタ補正逆投影法、逐次近似法等の単一波長用画像再構成アルゴリズムを用いて逆投影して画像(CT画像)を再構成する。
【0003】
上述のX線CT画像再構成方法において、撮影対象物にX線吸収係数が極端に大きい物質領域たとえば金属が存在すると、その物質領域を波長分布を有する照射X線が透過することによりその線質が変化するビームハードニング現象により、画像アーチファクトと呼ばれる放射状ノイズパターンがCT画像に発生する。この画像アーチファクトは医療現場、非破壊検査現場等で再構成画像の判別の妨げとなる。
【0004】
次に、ビームハードニング現象について、図40、図41、図42、図43を参照して説明する。
【0005】
図40は一般的なX線CT画像再構成装置を示す図である。
【0006】
図40において、1はX線発生器、2はX線検出器、3は撮影対象物、4は、X線発生器1、X線検出器2、撮影対象物3を制御する制御回路である。制御回路4はコンピュータで構成され、CPU、ROM、RAM等を含む。
【0007】
図41は図40の撮影対象物3の物質のX線吸収係数μを示す。
【0008】
図41に示すように、X線吸収係数μの分布が物質毎に異なるだけでなく、同一物質でもX線波長によってX線吸収係数μが異なる、つまり、エネルギー依存性を有する。
【0009】
図42は図40のCT画像再構成装置のX線エネルギースペクトルを示し、(A)はX線発生器1の発生X線の連続エネルギースペクトル(以下、単に、連続スペクトル)を示し、(B)は大きさL、X線吸収係数μ(E)の撮影対象物3を透過後のX線検出器2の検出X線の連続スペクトルを示す。
【0010】
図42に示すように、撮影対象物3の透過前後でX線の連続スペクトルは異なる。
【0011】
ところで、図40のX線CT画像再構成装置におけるCT画像再構成演算においては、X線吸収量のエネルギー依存性を考慮していない。つまり、μ一定の単一エネルギースペクトル(以下単に、単一スペクトル)のX線エネルギー照射つまり単一波長照射を前提としている。この理由は、X線連続スペクトルは装置、照射するX線等の種類、撮影条件等に依存するので、これに対応することが困難であるためである。しかし、実際のCT撮影は通常のX線のために波長分布を有する。この結果、単一波長用CT画像再構成時の演算によって矛盾が発生し、図43に示すような画像アーチファクトが発生する。
【0012】
画像アーチファクトを低減する第1の従来のX線CT画像再生方法は投影データのシノグラムから高X線吸収係数の物質領域たとえば金属領域を抽出して仮想投影によりその物質領域の投影データを改変して修正投影データを演算し、この修正投影データを逆投影して再構成画像(CT画像)を再構成している(参照:特許文献1)。
【0013】
上述の第1の従来のX線CT画像構成方法においては、厳しい画像アーチファクトを低減することができるが、軟組織部分の抽出が難しく投影データの修正が不充分であった。
【0014】
画像アーチファクトを低減する第2の従来のX線CT画像構成方法は、投影条件を変更して複数回CT撮影を行って複数の投影データを得てX線連続スペクトルの種類を限定し、X線吸収量を求めている(参照:特許文献2、非特許文献1)。
【先行技術文献】
【特許文献】
【0015】
【特許文献1】特開2006−167161号公報
【特許文献2】特開2008−241376号公報
【非特許文献】
【0016】
【非特許文献1】李根旭、小関道彦、木村仁、伊能教夫:X線CT画像の画質向上に関する研究(X線スペクトルを考慮した画像再構成手法の提案)、計測自動制御学会 第24回センシングフォーラム資料、pp.107-111,2007
【発明の概要】
【発明が解決しようとする課題】
【0017】
しかしながら、上述の第2の従来のX線CT画像再構成方法においては、複数のX線連続スペクトルを必要とするので、CT撮影回数が増加すると共に、撮影条件を変更することによる異なるX線連続スペクトルを得ることが難しいという課題がある。
【課題を解決するための手段】
【0018】
上述の課題を解決するために本発明は、X線CT画像再構成方法において、発生X線連続エネルギースペクトルを既知情報として記憶し、撮影対象物の単一エネルギースペクトルのX線投影データを生成し、補正量を初期化し、単一エネルギースペクトルのX線投影データを補正量で補正して撮影対象物に含まれる少なくとも2つの物質のX線吸収係数のエネルギー依存性を考慮して撮影対象物に単一エネルギースペクトルのX線を照射した場合に相当する補正X線投影データを得、補正X線投影データを単一波長用画像再構成アルゴリズムを用いて逆投影演算してCT画像を再構成し、各物質を区別するための少なくとも1つのしきい値を用いてCT画像から各物質の領域を抽出してこれらの領域のX線経路の長さを演算し、発生X線連続エネルギースペクトルのエネルギー、各物質の領域の長さ及び各物質のX線吸収係数を用いて連続エネルギースペクトルのX線投影データを演算し、単一エネルギースペクトルのX線投影データと連続エネルギースペクトルのX線投影データとの差により補正量を演算するものである。すなわち、従来の単一波長用の画像再構成アルゴリズムを用いて得られるCT画像(元CT画像)に対して、推定される撮影対象物に対して単一エネルギースペクトルのX線を照射したときX線投影データを推定する。この推定には、元CT画像に少なくとも1つのしきい値を用いて少なくとも2つの物質の領域を抽出し、各X線経路の長さを演算する。次いで、既知情報のX線連続エネルギースペクトルのエネルギー、物質の領域の長さ及びX線吸収係数を用いて実際に照射された連続エネルギースペクトルのX線投影データを補正演算し、単一エネルギースペクトルのX線投影データに近い値を得るようにしたものである。このようにして、少なくとも1つのしきい値で区分された少なくとも2つの物質が及ぼすビームハードニングの度合を個々に換算して補正量で表わす。
【0019】
また、上述の補正から補正量演算を繰返して補正量が所定範囲内に収束させる。この場合、しきい値が未知であれば、CT画像のヒストグラムによりしきい値を推定してしきい値の数を増加させる。このようにして、補正量が更新され、画像アーチファクトが徐々に減少することで物質の領域のX線経路の長さが正確となる。
【発明の効果】
【0020】
本発明によれば、単一エネルギースペクトルのX線投影データは変更されないので、CT撮影は1回で済み、しかも、補正を繰返すことによって画像アーチファクトを低減することができる。
【図面の簡単な説明】
【0021】
【図1】本発明に係るX線CT画像再構成方法の第1の実施の形態を説明するためのフローチャートである。
【図2】図1のステップ101のX線連続スペクトルの例を示すグラフである。
【図3】図1のステップ106のしきい値を説明する図であって、(A)はCT画像、(B)は画素値ヒストグラムである。
【図4】図1のステップ108の補正量を説明するためのグラフである。
【図5】X線発生器の発生X線の連続スペクトルを示す表である。
【図6】生体骨のX線吸収係数を示す表である。
【図7】図1の第1のシミュレーション例に用いられる撮影対象物としての骨円柱モデルを示す図である。
【図8】図7の補正前の骨円柱モデルを示し、(A)はCT画像、(B)は(A)センタライン上のCT値を示すグラフである。
【図9】図8の(A)のCT画像から物質の領域抽出を示すCT画像である。
【図10】図7の1回目補正の骨円柱モデルを示し、(A)はCT画像、(B)は(A)のセンタライン上のCT値を示すグラフである。
【図11】鉄のX線吸収係数を示す表である。
【図12】アクリル樹脂のX線吸収係数を示す表である。
【図13】図1の第2のシミュレーション例に用いられる撮影対象物としての複数同種ピンモデルを示す図であって、(A)は斜視像、(B)は上面像である。
【図14】図13の補正前の複数同種ピンモデルを示し、(A)はCT画像、(B)は(A)センタライン上のCT値を示すグラフである。
【図15】図14の(A)のCT画像から物質の領域抽出を示し、(A)は鉄の領域抽出を示すCT画像、(B)はアクリル樹脂の領域抽出を示すCT画像である。
【図16】図13の1回目補正の複数同種ピンモデルを示し、(A)はCT画像、(B)は(A)のセンタライン上のCT値を示すグラフである。
【図17】アルミニウムのX線吸収係数を示す表である。
【図18】図1の第3のシミュレーション例に用いられる撮影対象物としての複数異種ピンモデルを示す上面像である。
【図19】図18の補正前の複数異種ピンモデルを示し、(A)はCT画像、(B)は(A)センタライン上のCT値を示すグラフである。
【図20】図19の画素値ヒストグラムを示すグラフである。
【図21】図19の(A)のCT画像から物質の領域抽出を示し、(A)は鉄の領域抽出を示すCT画像、(B)はアルミニウムの領域抽出を示すCT画像、(C)はアクリル樹脂の領域抽出を示すCT画像である。
【図22】図18の1回目補正の複数異種ピンモデルを示し、(A)はCT画像、(B)は(A)の画素値ヒストグラムを示すグラフである。
【図23】図18の2回目補正の複数異種ピンモデルを示し、(A)はCT画像、(B)は(A)の画素値ヒストグラムを示すグラフである。
【図24】図18の3回目補正の複数異種ピンモデルを示し、(A)はCT画像、(B)は(A)の画素値ヒストグラムを示すグラフである。
【図25】図1の第3のシミュレーションの補正量の収束を示すグラフである。
【図26】図1の第1の実験例を示し、(A)は補正前のCT画像、(B)は補正後のCT画像である。
【図27】図1の第2の実験例を示し、(A)は補正前のCT画像、(B)は補正後のCT画像である。
【図28】図1の第3の実験例を示し、(A)は鉄釘―米粒モデルの斜視像、(B)は撮影時のX線スペクトルを示すグラフである。
【図29】図1の第3の実験例を示し、(A)は補正前のCT画像、(B)は補正後のCT画像である。
【図30】第3の実験例において鉄のみ抽出による補正後のCT画像である。
【図31】本発明に係るX線CT画像再構成方法の第2の実施の形態を説明するためのフローチャートである。
【図32】図31のステップ202、203における物質鉄、アルミニウム、アクリル樹脂のX線吸収係数の例を示すグラフである。
【図33】図31の第1のシミュレーション例の1回目補正後のCT画像である。
【図34】図31の第2のシミュレーション例を示し、(A)は1回目補正後のCT画像、(B)は補正後のCT画像である。
【図35】本発明に係るX線CT画像再構成方法の第3の実施の形態を説明するためのフローチャートである。
【図36】図35の第1のシミュレーション例を示し、(A)、(B)、(C)、(D)、(E)、(F)は、補正前、1回目補正後、2回目補正後、3回目補正後、4回目補正後、5回目補正後のCT画像である。
【図37】図35の第2のシミュレーション例を示し、(A)、(B)、(C)、(D)、(E)は、補正前、1回目補正後、2回目補正後、3回目補正後、4回目補正後のCT画像である。
【図38】図35の第1のシミュレーション例に用いられる撮影対象物としての生体組織モデルを示す図である。
【図39】図38のCT画像を示し、(A)は補正前のCT画像、(B)は補正後のCT画像である。
【図40】一般的なX線CT画像再構成装置を示す図である。
【図41】図40の撮影対象物の物質のX線吸収係数を示すグラフである。
【図42】図40のX線連続スペクトルを示すグラフである。
【図43】画像アーチファクトの例を示すCT画像である。
【発明を実施するための形態】
【0022】
図1は本発明に係るX線CT画像再構成方法の第1の実施の形態を説明するためのフローチャートである。図1のフローチャートは図40のX線CT画像再構成装置によって実行される。
【0023】
始めに、ステップ101では、X線発生器1のX線連続スペクトルを予め記憶する。このX線発生器1のX線管の種類、撮影条件によって定まるX線連続スペクトルは専用の機器で予め計測できる。尚、X線連続スペクトルは実際には離散値としてEi及びその光子数n(Ei)で記憶される。
【0024】
次に、ステップ102では、後述の補正量ht及び補正回数tを初期化する。つまり、
h0←0
t←0
ここで、h0は投影データの数だけ存在するので、h0はその総称を表している。
【0025】
次にステップ103では、X線検出器2の検出信号により単一スペクトルつまり波長iでの単一スペクトル投影データ(X線吸収量)gを生成する。尚、波長iはX線管電圧から計算される代表値を採用する。
【0026】
ステップ103におけるX線吸収量gを理論的に考察すると、X線吸収量gは数1の式で表される。
【数1】
但し、Iin:物質kの入射時のエネルギー
Iout:物質kの透過後のエネルギー
Lk:物質kのX線経路の長さ
μik:波長i及び物質kに依存するX線吸収係数
X線が透過する物質がn種類の物質で構成されていれば、数1の式は数2の式となる
【数2】
数2の式を対数で表すと、数3の式となる。
【数3】
【0027】
次に、ステップ104では、単一スペクトル投影データgに補正量htを減算して補正投影データg’を求める。
g’= g - ht
尚、ここでも、g’は投影データの数だけ存在するので、その総称を表している。
【0028】
次に、ステップ105では、補正投影データg’をフィルタ補正逆投影法等により逆投影してCT画像を再構成する。この場合の逆投影法は、従来の単一スペクトルのX線投影データを再構成するための単一波長用画像再構成アルゴリズムを用いる。
【0029】
次に、ステップ106では、既知のしきい値TH1、TH2、…を用いてCT画像から物質kのX線経路の領域の長さLkを演算する。たとえば、図3の(A)に示すCT画像が図3の(B)に示すごとく3種類の物質、つまり、鉄、アルミニウム、アクリル樹脂を含んでいれば、ステップ106において、TH1を画素値14000、TH2を画素値600に予め設定しておく。尚、対象撮影物3が複数の未知の物質を含んでいる場合には、ステップ106において、図3の(B)の画素値ヒストグラムを演算し、画素値の集中箇所、たとえば、図3の(B)における画素値-900、1800、8000、20000の集中箇所を判別し、その間にしきい値を推定して設定する。尚、このしきい値推定アルゴリズムは容易に得ることができる。
【0030】
次に、ステップ107では、各単一スペクトル投影データgに対応する連続スペクトル投影データfを演算する。
【0031】
連続スペクトル投影データfはステップ101にて記憶されたX線発生器1のX線連続スペクトルのエネルギーEi、物質kのX線吸収係数μik及びステップ106にて演算された物質kのX線経路の長さLkに基づいて演算される。つまり、物質kの入射時のX線エネルギーIinは数4の式で表わされる。
【数4】
また、物質kの透過後のX線エネルギーIoutは数5の式で表わされる。
【数5】
数4、数5の式から数6の式が得られる。
【数6】
数6の式を対数で表わすと、連続スペクトル投影データ(X線吸収量)fは数7の式で表わされる。
【数7】
【0032】
次に、ステップ108では、
t ← t + 1
とする。
【0033】
次に、ステップ109では、t回目の補正量htを演算する。すなわち、
ht = g - f
【0034】
すなわち、図4に示すように、物質kのX線経路の長さLkが同一であっても、単一スペクトル投影データ(X線吸収量)gと連続スペクトル投影データ(X線吸収量)fとでは値が補正量htだけ異なっている。この補正量htだけ異なっている分、この補正量htが画像アーチファクトとして現われている。ステップ106における画像アーチファクトが存在する投影データg’からの物質kのX線経路の長さLkにも僅差が存在する。従って、ステップ104において単一スペクトル投影データgをこの補正量htで減算して補正することにより連続スペクトル投影データfに近づけることができ、長さLkも安定的収束し、これに伴い、画像アーチファクト低減効果が生ずる。
【0035】
ステップ110はステップ104〜109のフローを補正量htが所定範囲に収束して|ht-ht-1|がε未満になるまで繰返す。|ht-ht-1|がε未満となったときにステップ111にてこのルーチンは終了する。尚、ここでも、|ht-ht-1|は投影データの数だけ存在し、その総和を意味する。
【0036】
図1のフローチャートで重要なのは、1回のCT撮影で得られた単一スペクトル投影データgは変更されず、その補正量htのみを更新し、この更新によって画像アーチファクトを徐々に低減することで物質kのX線経路の長さLkの演算も正確になっていくことである。
【0037】
図1のフローチャートによるシミュレーション例を説明する。つまり、この場合、制御回路(コンピュータ)4上でCT撮影する仮想投影手法を用い、X線発生器1のX線連続スペクトルは図5に示す9個の離散値で近似する。
【0038】
第1のシミュレーション例は図6に示すX線吸収係数μikを有する物質つまり生体骨を用いた図7に示す骨円柱モデルに適用する。つまり、物質k及びそのX線吸収係数μikは既知とする。この場合、骨円柱モデルは均一であり、かつ形状が単純であるので、ビームハードニングの影響が小さく、従って、画像アーチファクトが少ないと考えられるケースである。
【0039】
第1のシミュレーション例においては、ステップ101において、図5のX線連続スペクトル及び図6のX線吸収係数μikを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。この場合、図8の(A)に示すCT画像が得られ、一見すると画像アーチファクトは出現していないが、図8の(B)に示すように、CT画像のセンタライン上のCT値をプロットすると、円柱中央部の画素で輝度が上昇し、円柱周辺部の画素で輝度が低下していることが分かる。つまり、円柱周辺部の画素のCT値が真値1240より大きくなっている。このように、単純な形状でもビームハードニングの影響で適切なX線吸収量が得られず、CT値の定量性が喪失してしまう。
【0040】
次に、図8の(A)のCT画像に1回目補正を行う。つまり、ステップ106において、図8の(A)のCT画像の画素に対して既知のしきい値TH1たとえば500を用いて図9に示す物質生体骨の領域を抽出する。図9に示すごとく、物質生体骨の領域は正確に抽出されている。この抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ107において、その長さLk、図5のX線連続スペクトルのエネルギーEi、図6のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データから図10の(A)に示すCT画像を再構成する。
【0041】
図10の(A)の1回目補正のCT画像についてセンタライン上の画素をプロットすると、図10の(B)に示すごとく、円柱周辺部の画素の輝度が上昇して円柱内部で均一となり、50keVのCT真値1240に近づいていることが分かる。
【0042】
第2のシミュレーション例は図11、図12に示すX線吸収係数μikを有する物質つまり鉄、アクリル樹脂を用いた図13に示す複数同種ピンモデルに適用する。類似モデルとしては、電子部品、半導体の非破壊検査の際に、金属部品が樹脂中に埋め込められていると、金属部品の周囲に画像アーチファクトが発生して金属部品の周囲の樹脂部分の形状、組成が判別できなくなる場合がある。また、人体にインプラント金属が埋め込められていると、インプラント金属の周囲に画像アーチファクトが発生してインプラント金属の周囲の筋肉部分の形状、組成が判別できなくなる場合がある。尚、この場合も、物質k及びそのX線吸収係数μikは既知とする。この複数同種ピンモデルは形状が少し複雑であるので、ビームハードニングの影響が大きく、従って、画像アーチファクトがある程度大きいと考えられるケースである。
【0043】
第2のシミュレーション例においては、ステップ101において、図5のX線連続スペクトルのエネルギーEi及び図11、図12のX線吸収係数μikを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。この場合、図14の(A)に示すCT画像が得られ、強い画像アーチファクトが出現している。図14の(B)に示すように、CT画像のセンタライン上のCT値をプロットすると、ピン中央部の画素で輝度が急上昇し、ピン周辺部の画素で輝度が急低下していることが分かる。つまり、ピン周辺部の画素のCT値が真値20000より大きくなっている。このように、複雑な形状では強いビームハードニングの影響で適切なX線吸収量が得られず、CT値の定量性が喪失してしまう。
【0044】
次に、図14の(A)のCT画像に1回目補正を行う。つまり、ステップ106において、図14の(A)のCT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて図15に示す物質鉄及びアクリル樹脂の各領域を抽出する。この場合、図15の(A)に示すごとく、物質鉄の領域は正確に抽出されているが、図15の(B)に示すごとく、物質アクリル樹脂の領域は物質鉄に起因するアーチファクトの影響で不正確である。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ107において、これらの長さLk、図5のX線連続スペクトルのエネルギーEi、図11、図12のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図16の(A)に示すCT画像を再構成する。
【0045】
図16の(A)の1回目補正のCT画像についてセンタライン上画素をプロットすると、図16の(B)に示すごとく、ピン周辺部の画素の輝度が上昇してピン内部で均一となり、50keVのCT真値20000に近づいていることが分かる。また、同時に、物質アクリル樹脂の領域の画素の輝度も均一となり、50keVのCT真値-680に近づいていることが分かる。
【0046】
第3のシミュレーション例は図11、図12、図17に示すX線吸収係数μikを有する物質つまり鉄、アクリル樹脂、アルミニウムを用いた図18に示す複数異種ピンモデルに適用する。尚、この場合も、物質k及びそのX線吸収係数μikは既知とする。この複数異種ピンモデルは複数の高X線吸収係数を有する物質を有するので、ビームハードニングの影響が非常に大きく、従って、画像アーチファクトが大きいと考えられるケースである。
【0047】
第3のシミュレーション例においては、ステップ101において、図5のX線連続スペクトルのエネルギーEi及び図11、図12、図17のX線吸収係数μikを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。この場合、図19の(A)に示すCT画像が得られ、物質鉄に起因する強い画像アーチファクトが出現している。図19の(B)に示すように、CT画像のセンタライン上のCT値をプロットすると、ピン中央部の画素で輝度が急上昇し、ピン周辺部の画素で輝度が急低下していることが分かる。つまり、ピン周辺部の画素のCT値が真値20000あるいは1800より大きくなっている。このときの画素ヒストグラムを図20に示すと、物質鉄の分布は真値20000から離れた23000程度の箇所にある。従って、図1のステップ106におけるしきい値TH1を8000と推定して設定し、物質鉄の領域を抽出できる。尚、この場合、物質鉄が未知であっても、1つのしきい値TH1が推定され設定される。この結果、図21の(A)に示すごとく、物質鉄の領域を正確に抽出できるが、図21の(B)、(C)に示すごとく、物質アルミニウム、アクリル樹脂の各領域は物質鉄に起因するアーチファクトの影響で正確には抽出できない。
【0048】
次に、図19の(A)のCT画像に1回目補正を行う。つまり、上述の抽出領域のX線経路毎に長さLkを演算する。次いで、これらの長さLk、図5のX線連続スペクトルのエネルギーEi、図11、図12、図17のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図22の(A)に示すCT画像を再構成する。図22の(A)のCT画像においては、画像アーチファクトが残り、補正は不充分である。この理由は、物質鉄に起因する画像アーチファクトの画素値が物質アルミニウムの画素値とほぼ同一であるために、画像アーチファクトの形状が物質アルミニウムが存在するものとして再構成されたものと考えられる。このときの画素ヒストグラムを図22の(B)に示すと、物質アルミニウム及びアクリル樹脂のピーク値がこれらの真値-680及び1800に現われている。従って、図1のステップ106におけるしきい値TH2を-680と1800との間たとえば1200と推定して設定すれば、物質アルミニウム、アクリル樹脂の各領域は正確に抽出できる。尚、この場合も、物質アルミニウム、アクリル樹脂が未知であっても、しきい値TH2が推定され設定される。従って、撮影対象物に含まれるn個の物質が未知である場合には、図1のステップ104〜110を繰返すことによりn−1個のしきい値が推定され設定されることになる。この結果、未知の物質の領域も正確に抽出できる。
【0049】
次に、図22の(A)のCT画像に2回目補正を行う。つまり、ステップ106において、図22の(A)のCT画像の画素に対して推定のしきい値TH1たとえば8000及びしきい値TH2たとえば1200を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。次いで、ステップ107において、これらの長さLk、図5のX線連続スペクトルのエネルギーEi、図11、図12、図17のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、2回目の補正量h2を演算し、補正量h2で単一スペクトル投影データgを補正し、補正された投影データg’から図23の(A)に示すCT画像を再構成する。
【0050】
次に、図23の(A)のCT画像に3回目補正を行う。つまり、ステップ106において、図23の(A)のCT画像の画素に対して推定のしきい値TH1たとえば8000及びしきい値TH2たとえば1200を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。次いで、ステップ107において、これらの長さLk、図5のX線連続スペクトルのエネルギーEi、図11、図12、図17のX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、3回目の補正量h3を演算し、補正量h3で単一スペクトル投影データgを補正し、補正された投影データg’から図23の(A)に示すCT画像を再構成する。
【0051】
図24の(A)の3回目補正のCT画像についてセンタライン上画素をプロットすると、図24の(B)に示すごとく、ピン周辺部の画素の輝度が上昇してピン内部で均一となり、50keVのCT真値20000もしくは1800に近づいていることが分かる。また、同時に、物質アクリル樹脂の領域の画素の輝度も均一となり、50keVのCT真値-680に近づいていることが分かる。
【0052】
第3のシミュレーション例においては、CT画像をたとえ4回目補正、5回目補正、…を行ってもCT画像はほとんど変わらない。つまり、補正回数が増加するにつれて、画像アーチファクトが低減し、物質鉄、アルミニウム、アクリル樹脂の領域がほぼ正確に抽出できる。従って、図25に示すごとく、補正回数が3回を超えて総和|ht-ht-1|がε以下となると、各物質がその真値に近づいており、CT画像の誤差は低減し、補正の効果が向上していることが確認できる。このとき、補正量と前回補正量との差|ht-ht-1|も0に近づき、補正量が一定値に収束していることが分かる。
【0053】
次に、図1のフローチャートによる実験例を説明する。
【0054】
図26は上述の第2のシミュレーション例に対応して行った第1の実験例である。撮影時の発生X線の連続スペクトルは不明であったので、図5の連続スペクトルの離散値で近似させた。また、物質鉄、アクリル樹脂のX線吸収係数μikは図11、図12に示す値を既知として用いた。さらに、物質鉄、アクリル樹脂の各領域の抽出のためのしきい値に適当な値を既知として用いた。この結果、補正前のCT画像は図26の(A)に示すごとくであったが、図1のフローチャートに基づく補正を行った結果、図26の(B)に示すCT画像が得られた。つまり、画像アーチファクトは大幅に減少し、補正の効果が認められた。
【0055】
図27は鉄の十字架形状モデルに対して行った第2の実験例である。この場合も、撮影時の発生X線の連続スペクトルは不明であったので、図5の連続スペクトルの離散値で近似させた。また、物質鉄、アクリル樹脂のX線吸収係数μikは図11、図12に示す値を既知として用いた。さらに、物質鉄、アクリル樹脂の各領域の抽出のためのしきい値に適当な値を既知として用いた。この結果、補正前のCT画像は図27の(A)に示すごとくであったが、図1のフローチャートに基づく補正を行った結果、図27の(B)に示すCT画像が得られた。やはり、画像アーチファクトは大幅に減少し、補正の効果が認められた。
【0056】
図28の(A)に示す9本の鉄釘を米粒の山へ埋込んだモデルに対して行ったものである。撮影時の管電圧130kVの発生X線の連続スペクトルは、図28の(B)の連続スペクトルの離散値で近似させた。また、物質鉄、米粒のX線吸収係数μikは図11、図12に示す値を既知として用いた。つまり、米粒はアクリル樹脂とみなした。さらに、物質鉄、米粒の各領域の抽出のためのしきい値に13500及び750を既知として用いた。この結果、補正前のCT画像は図29の(A)に示すごとくであった。つまり、鉄釘と鉄釘とを結ぶラインで強い画像アーチファクトが発生している。これに対し、図1のフローチャートに基づく補正を行った結果、図29の(B)に示すCT画像が得られた。やはり、画像アーチファクトは大幅に減少し、補正の効果が認められた。
【0057】
尚、第3の実験例において、しきい値13500だけ設定して鉄釘(鉄)のみの領域を抽出した場合には、図30に示すごとく、米粒が正確に抽出できなかった。この理由は、連続スペクトルにより、鉄釘(鉄)だけでなく、米粒部分でもX線吸収量の過吸収が発生しているものと考えられるからである。
【0058】
このように、実験例によれば、本発明はCT画像を用いた個体別モデリング、リバースエンジニアリング非破壊検査で有効であることが分かる。
【0059】
図31は本発明に係るX線CT画像再構成方法の第2の実施の形態を説明するためのフローチャートであって、対象撮影物3を構成する物質kの数が少なくかつこれら物質kのX線吸収係数μikが未知の場合に適用される。
【0060】
図31においては、図1において、ステップ101をステップ201に置換し、ステップ202、203を追加したものである。すなわち、ステップ201では、図5のX線連続スペクトルのエネルギーEiのみを記憶し、ステップ202、203において、X線連続スペクトルのエネルギーEi及び物質kのX線経路の長さLkに基づいて物質kのX線吸収係数μikを演算する。
【0061】
ステップ202、203の根拠を説明する。X線検出器2のある位置の投影データDは数8の式で表わされる。
【数8】
ここで、X線連続スペクトルのエネルギーEiが既知であるので、数8の式の投影データはX線吸収係数μik及びステップ106で定まる長さLkの関数となる。他方、X線吸収係数μikは数9の式で表わされる。
【数9】
数9の式の妥当性は、図32に示すように、物質鉄、アルミニウム、アクリル樹脂に対してak=134000、7000、200としたときに、数9の式により求められたX線吸収係数μikは実際のX線吸収係数μikにほぼ一致することから確認できる。
数9の式で数8の式を表わすと、数10の式となる。
【数10】
数10の式において、物質kの種類が未知数n個のn次方程式としてakを解くことができる。実際に検出されるデータはX線検出器2の列数たとえば512×撮影の回転方向数たとえば800個分入手できるので、物質kの種類数nはどんなに多くとも解くことができる。つまり、複数の方向から照射したX線ビーム毎に投影データD’が得られるので、数10の式がその数だけ成立する。この場合、数10の式の数はn以上であればakを求めることができる。但し、領域の長さLkは画像アーチファクトによる誤差を有するので、akは最小2乗法で求める。尚、このとき、求められたakに基づき、物質kを推定できる。
【0062】
従って、ステップ202では、数10の式のakを最小2乗法により演算する。次いで、ステップ203において、物質kのX線吸収係数μikを数9の式により演算する。
【0063】
図31のフローチャートによるシミュレーション例を説明する。つまり、この場合、制御回路(コンピュータ)4上でCT撮影する仮想投影手法を用い、X線発生器1のX線連続スペクトルは図5に示す9個の離散値で近似する。
【0064】
第1のシミュレーション例は2つの物質つまり鉄、アクリル樹脂を用いた図13に示す複数同種ピンモデルに適用する。但し、物質鉄、アクリル樹脂のX線吸収係数μikは未知とする。
【0065】
第1のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。
【0066】
次に、CT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202において、数10の式のakを最小2乗法により演算する。その結果、
a1 = 124651
a2 = 41
が得られた。尚、a(鉄)、a(アクリル樹脂)の真値は130000、200であるので、上述の結果値は真値に近いことが分かる。従って、a1 = 124651は鉄に相当、a2= 41はアクリル樹脂に相当すると推定する。次いで、ステップ203において、物質鉄及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図33に示すCT画像を再構成する。このように、X線吸収係数μikが未知である場合にも、画像アーチファクトが大幅に減少する。
【0067】
第2のシミュレーション例は3つの物質つまり鉄、アルミニウム及びアクリル樹脂を用いた図18に示す複数異種ピンモデルに適用する。但し、物質鉄、アルミニウム及びアクリル樹脂のX線吸収係数μikは未知とする。
【0068】
第2のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgからCT画像を再構成する。
【0069】
次に、CT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202において、数10の式のakを最小2乗法により演算する。その結果、
a1 = 127682
a2 = -1702
a3 = 416
が得られた。この場合、画像アーチファクトの影響によりアルミニウム相当の領域抽出が不適切であるために、a(アルミニウム)の真値7000から離れてしまった。また、X線吸収係数は負とならないので、
a1 = 127682
a2 = 0
a3 = 416
とする。次いで、ステップ203において、物質鉄、アルミニウム及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図34の(A)に示すCT画像を再構成する。このように、1回目補正により、不充分であるが、画像アーチファクトが減少する。
【0070】
次に、図34の(A)のCT画像に2回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄、アルミニウム及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202において、数10の式のakを最小2乗法により演算する。その結果、
a1 = 123476
a2 = 6962
a3 = 21
が得られた。尚、a(鉄)、a(アルミニウム)、a(アクリル樹脂)の真値は130000、7000、200であるので、上述の結果値は真値に近いことが分かる。従って、a1 = 123476は鉄、a2= 6962はアルミニウム、a3 = 21はアクリル樹脂に相当すると推定する。次いで、ステップ203において、物質鉄、アルミニウム及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図34の(B)に示すCT画像を再構成する。このように、2回目補正により、画像アーチファクトが大幅に減少する。
【0071】
図35は本発明に係るX線CT画像再構成方法の第3の実施の形態を説明するためのフローチャートであって、対象撮影物3を構成する物質kの数が多くかつこれら物質kのX線吸収係数μikが未知の場合に適用される。
【0072】
図35においては、図31において、ステップ202をステップ202Aに置換したものである。すなわち、ステップ201Aでは、図31のフローチャートにおけるakを
ak = 9.873・THn + 9873
= 9.873・(THn+ 1000)
但し、THn(n=1,2,…)は暫定的なしきい値、
により演算する。ここで、値9.873、9873は、空気のCT値が-1000である場合にa(空気)=0、鉄のCT値12572の場合にa(鉄)=134000となるように線形に規定したものである。ここで、鉄のCT値12572は実際のCT画像の最大CT値である。但し、9.873、9873は与えられた2点が線形となるように、他の値αk、βkとなし得る。この場合も、求められたakに基づき、物質kを推定できる。
【0073】
図35のフローチャートによるシミュレーション例を説明する。つまり、この場合、制御回路(コンピュータ)4上でCT撮影する仮想投影手法を用い、X線発生器1のX線連続スペクトルは図5に示す9個の離散値で近似する。
【0074】
第1のシミュレーション例は2つの物質つまり鉄、アクリル樹脂を用いた図13に示す複数同種ピンモデルに適用する。但し、物質鉄、アクリル樹脂のX線吸収係数μikは未知とする。
【0075】
第1のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgから図36の(A)に示すCT画像を再構成する。
【0076】
次に、CT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202Aにおいて、akを演算する。その結果、
a(鉄相当)= 128349
a(アクリル樹脂相当)= 8885
が得られた。次いで、ステップ203において、物質鉄及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図36の(B)に示すCT画像を再構成する。同様にして、2回目補正、3回目補正、4回目補正、5回目補正を行うと、図36の(C)、(D)、(E)、(F)に示すCT画像が得られる。計算回数は増加するが、画像アーチファクトが徐々に減少する。
【0077】
第2のシミュレーション例は3つの物質つまり鉄、アルミニウム、アクリル樹脂を用いた図18に示す複数異種ピンモデルに適用する。但し、物質鉄、アルミニウム、アクリル樹脂のX線吸収係数μikは未知とする。
【0078】
第2のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgから図37の(A)に示すCT画像を再構成する。
【0079】
次に、CT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202Aにおいて、akを演算する。その結果、
a(鉄相当)= 128349
a(アルミニウム相当)= 8885
a(アクリル樹脂相当)= 8885
が得られた。次いで、ステップ203において、物質鉄及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’から図37の(B)に示すCT画像を再構成する。同様にして、2回目補正、3回目補正、4回目補正、5回目補正を行うと、図37の(C)、(D)、(E)に示すCT画像が得られる。計算回数は増加するが、やはり、画像アーチファクトが徐々に減少する。
【0080】
第3のシミュレーション例は図38に示す生体組織モデルに適用する。但し、生体組織のX線吸収係数μikは未知とする。
【0081】
第3のシミュレーション例においては、ステップ201において、図5のX線連続スペクトルのエネルギーEiを記憶する。そして、ステップ103において、これらに基づいて単一スペクトル投影データgを生成し、ステップ105において、単一スペクトル投影データgから図39の(A)に示すCT画像を再構成する。
【0082】
次に、図39の(A)に示すCT画像に1回目補正を行う。つまり、ステップ106において、CT画像の画素に対して既知のしきい値TH1たとえば1300及びしきい値TH2たとえば900を用いて物質鉄及びアクリル樹脂の各領域を抽出する。これらの抽出領域のX線経路毎に長さLkを演算する。次いで、ステップ202Aにおいて、akを演算する。その結果、
a(鉄相当)= 128349
a(アクリル樹脂相当)= 8885
が得られた。次いで、ステップ203において、物質鉄及びアクリル樹脂相当のX線吸収係数μikを数9の式を用いて演算する。次いで、ステップ107において、上述の長さLk、図5のX線連続スペクトルのエネルギーEi、ステップ203にて得られたX線吸収係数μikを用いて連続スペクトル投影データfを演算し、次いで、ステップ109において、1回目の補正量h1を演算し、補正量h1で単一スペクトル投影データgを補正し、補正された投影データg’からCT画像を再構成する。同様にして、複数回補正を行うと、図39の(B)に示すCT画像が得られる。計算回数は増加するが、やはり、画像アーチファクトが徐々に減少する。
【符号の説明】
【0083】
1:X線発生器
2:X線検出器
3:対象撮影物
4:制御回路
Iin:入射時X線エネルギー
Iout:透過後X線エネルギー
μik:X線吸収係数
g:単一スペクトル投影データ
f:連続スペクトル投影データ
ht:補正量
【特許請求の範囲】
【請求項1】
発生X線連続エネルギースペクトルを既知情報として記憶する記憶段階と、
撮影対象物の単一エネルギースペクトルのX線投影データを生成するX線投影データ生成段階と、
補正量を初期化する初期化段階と、
前記単一エネルギースペクトルのX線投影データを前記補正量で補正して前記撮影対象物に含まれる少なくとも2つの物質のX線吸収係数のエネルギー依存性を考慮して前記撮影対象物に単一エネルギースペクトルのX線を照射した場合に相当する補正X線投影データを得る補正段階と、
前記補正X線投影データを単一波長用画像再構成アルゴリズムを用いて逆投影演算してCT画像を再構成するCT画像再構成段階と、
前記各物質を区別するための少なくとも1つのしきい値を用いて前記CT画像から前記各物質の領域を抽出して該領域のX線経路の長さを演算する長さ演算段階と、
前記発生X線連続エネルギースペクトルのエネルギー、前記各物質の領域の長さ及び前記各物質のX線吸収係数を用いて連続エネルギースペクトルのX線投影データを演算するX線投影データ演算段階と、
前記単一エネルギースペクトルのX線投影データと前記連続エネルギースペクトルのX線投影データとの差により前記補正量を演算する補正量演算段階と
を具備するX線CT画像再構成方法。
【請求項2】
前記物質が既知であり、さらに、前記補正段階、前記CT画像再構成段階、前記長さ演算段階、前記X線投影データ演算段階及び前記補正量演算段階を繰返して前記補正量が所定範囲に収束するようにする段階を具備する請求項1に記載のX線CT画像再構成方法。
【請求項3】
前記物質が未知であり、さらに、前記CT画像のヒストグラムから前記しきい値を推定するしきい値推定段階を具備する請求項1に記載のX線CT画像再構成方法。
【請求項4】
さらに、前記補正段階、前記CT画像再構成段階、前記しきい値推定段階、前記長さ演算段階、前記X線投影データ演算段階及び前記補正量演算段階を繰返して前記補正量が所定範囲に収束するようにする段階を具備する請求項3に記載のX線CT画像再構成方法。
【請求項5】
前記各物質のX線吸収係数は予め設定された請求項1に記載のX線CT画像再構成方法。
【請求項6】
さらに、
前記各物質kの係数akを
【数1】
を用いて演算する段階と、
該演算された係数ak及び前記発生X線連続エネルギースペクトルのエネルギーEiを用いて前記物質のX線吸収係数μikを
【数2】
により演算する段階と
を具備する請求項1に記載のX線CT画像再構成方法。
【請求項7】
さらに、
前記各物質kの係数akを
ak=αk・THn(n=1,2,…)+β
但し、THnは前記しきい値、αk、βkは定数
により演算する段階と、
該演算された係数ak及び前記発生X線連続エネルギースペクトルのエネルギーEiを用いて前記各物質kのX線吸収係数μikを
【数2】
により演算する段階と
を具備する請求項1に記載のX線CT画像再構成方法。
【請求項8】
さらに、前記係数akにより前記物質kを推定する段階を具備する請求項6あるいは7に記載のX線CT画像再構成方法。
【請求項9】
発生X線連続エネルギースペクトルを既知情報として記憶する記憶手順と、
撮影対象物の単一エネルギースペクトルのX線投影データを生成するX線投影データ生成手順と、
補正量を初期化する初期化手順と、
前記単一エネルギースペクトルのX線投影データを前記補正量で補正して前記撮影対象物に含まれる少なくとも2つの物質のX線吸収係数のエネルギー依存性を考慮して前記撮影対象物に単一エネルギースペクトルのX線を照射した場合に相当する補正X線投影データを得る補正手順と、
前記補正X線投影データを単一波長用画像再構成アルゴリズムを用いて逆投影演算してCT画像を再構成するCT画像再構成手順と、
前記各物質を区別するための少なくとも1つのしきい値を用いて前記CT画像から前記各物質の領域を抽出して該領域のX線経路の長さを演算する長さ演算手順と、
前記発生X線連続エネルギースペクトルのエネルギー、前記各物質の領域の長さ及び前記各物質のX線吸収係数を用いて連続エネルギースペクトルのX線投影データを演算するX線投影データ演算手順と、
前記単一エネルギースペクトルのX線投影データと前記連続エネルギースペクトルのX線投影データとの差により前記補正量を演算する補正量演算手順と
を具備するX線CT画像再構成プログラム。
【請求項10】
前記物質が既知であり、さらに、前記補正手順、前記CT画像再構成手順、前記長さ演算手順、前記X線投影データ演算手順及び前記補正量演算手順を繰返して前記補正量が所定範囲に収束するようにする手順を具備する請求項9に記載のX線CT画像再構成プログラム。
【請求項11】
前記物質が未知であり、さらに、前記CT画像のヒストグラムから前記しきい値を推定するしきい値推定手順を具備する請求項9に記載のX線CT画像再構成プログラム。
【請求項12】
さらに、前記補正手順、前記CT画像再構成手順、前記しきい値推定手順、前記長さ演算手順、前記X線投影データ演算手順及び前記補正量演算手順を繰返して前記補正量が所定範囲に収束するようにする手順を具備する請求項11に記載のX線CT画像再構成プログラム。
【請求項13】
前記各物質のX線吸収係数は予め設定された請求項9に記載のX線CT画像再構成プログラム。
【請求項14】
さらに、
前記物質kの係数akを
【数1】
を用いて演算する手順と、
該演算された係数ak及び前記発生X線連続エネルギースペクトルのエネルギーEiを用いて前記物質kのX線吸収係数μikを
【数2】
により演算する手順と
を具備する請求項9に記載のX線CT画像再構成プログラム。
【請求項15】
さらに、
前記物質kの係数akを
ak=αk・THn(n=1,2,…)+βk
但し、THnは前記しきい値、αk、βkは定数
により演算する手順と、
該演算された係数ak及び前記発生X線連続エネルギースペクトルのエネルギーEiを用いて前記各物質kのX線吸収係数μikを
【数2】
により演算する手順と
を具備する請求項7に記載のX線CT画像再構成プログラム。
【請求項16】
さらに、前記係数akにより前記物質kを推定する段階を具備する請求項14あるいは15に記載のX線CT画像再構成プログラム。
【請求項1】
発生X線連続エネルギースペクトルを既知情報として記憶する記憶段階と、
撮影対象物の単一エネルギースペクトルのX線投影データを生成するX線投影データ生成段階と、
補正量を初期化する初期化段階と、
前記単一エネルギースペクトルのX線投影データを前記補正量で補正して前記撮影対象物に含まれる少なくとも2つの物質のX線吸収係数のエネルギー依存性を考慮して前記撮影対象物に単一エネルギースペクトルのX線を照射した場合に相当する補正X線投影データを得る補正段階と、
前記補正X線投影データを単一波長用画像再構成アルゴリズムを用いて逆投影演算してCT画像を再構成するCT画像再構成段階と、
前記各物質を区別するための少なくとも1つのしきい値を用いて前記CT画像から前記各物質の領域を抽出して該領域のX線経路の長さを演算する長さ演算段階と、
前記発生X線連続エネルギースペクトルのエネルギー、前記各物質の領域の長さ及び前記各物質のX線吸収係数を用いて連続エネルギースペクトルのX線投影データを演算するX線投影データ演算段階と、
前記単一エネルギースペクトルのX線投影データと前記連続エネルギースペクトルのX線投影データとの差により前記補正量を演算する補正量演算段階と
を具備するX線CT画像再構成方法。
【請求項2】
前記物質が既知であり、さらに、前記補正段階、前記CT画像再構成段階、前記長さ演算段階、前記X線投影データ演算段階及び前記補正量演算段階を繰返して前記補正量が所定範囲に収束するようにする段階を具備する請求項1に記載のX線CT画像再構成方法。
【請求項3】
前記物質が未知であり、さらに、前記CT画像のヒストグラムから前記しきい値を推定するしきい値推定段階を具備する請求項1に記載のX線CT画像再構成方法。
【請求項4】
さらに、前記補正段階、前記CT画像再構成段階、前記しきい値推定段階、前記長さ演算段階、前記X線投影データ演算段階及び前記補正量演算段階を繰返して前記補正量が所定範囲に収束するようにする段階を具備する請求項3に記載のX線CT画像再構成方法。
【請求項5】
前記各物質のX線吸収係数は予め設定された請求項1に記載のX線CT画像再構成方法。
【請求項6】
さらに、
前記各物質kの係数akを
【数1】
を用いて演算する段階と、
該演算された係数ak及び前記発生X線連続エネルギースペクトルのエネルギーEiを用いて前記物質のX線吸収係数μikを
【数2】
により演算する段階と
を具備する請求項1に記載のX線CT画像再構成方法。
【請求項7】
さらに、
前記各物質kの係数akを
ak=αk・THn(n=1,2,…)+β
但し、THnは前記しきい値、αk、βkは定数
により演算する段階と、
該演算された係数ak及び前記発生X線連続エネルギースペクトルのエネルギーEiを用いて前記各物質kのX線吸収係数μikを
【数2】
により演算する段階と
を具備する請求項1に記載のX線CT画像再構成方法。
【請求項8】
さらに、前記係数akにより前記物質kを推定する段階を具備する請求項6あるいは7に記載のX線CT画像再構成方法。
【請求項9】
発生X線連続エネルギースペクトルを既知情報として記憶する記憶手順と、
撮影対象物の単一エネルギースペクトルのX線投影データを生成するX線投影データ生成手順と、
補正量を初期化する初期化手順と、
前記単一エネルギースペクトルのX線投影データを前記補正量で補正して前記撮影対象物に含まれる少なくとも2つの物質のX線吸収係数のエネルギー依存性を考慮して前記撮影対象物に単一エネルギースペクトルのX線を照射した場合に相当する補正X線投影データを得る補正手順と、
前記補正X線投影データを単一波長用画像再構成アルゴリズムを用いて逆投影演算してCT画像を再構成するCT画像再構成手順と、
前記各物質を区別するための少なくとも1つのしきい値を用いて前記CT画像から前記各物質の領域を抽出して該領域のX線経路の長さを演算する長さ演算手順と、
前記発生X線連続エネルギースペクトルのエネルギー、前記各物質の領域の長さ及び前記各物質のX線吸収係数を用いて連続エネルギースペクトルのX線投影データを演算するX線投影データ演算手順と、
前記単一エネルギースペクトルのX線投影データと前記連続エネルギースペクトルのX線投影データとの差により前記補正量を演算する補正量演算手順と
を具備するX線CT画像再構成プログラム。
【請求項10】
前記物質が既知であり、さらに、前記補正手順、前記CT画像再構成手順、前記長さ演算手順、前記X線投影データ演算手順及び前記補正量演算手順を繰返して前記補正量が所定範囲に収束するようにする手順を具備する請求項9に記載のX線CT画像再構成プログラム。
【請求項11】
前記物質が未知であり、さらに、前記CT画像のヒストグラムから前記しきい値を推定するしきい値推定手順を具備する請求項9に記載のX線CT画像再構成プログラム。
【請求項12】
さらに、前記補正手順、前記CT画像再構成手順、前記しきい値推定手順、前記長さ演算手順、前記X線投影データ演算手順及び前記補正量演算手順を繰返して前記補正量が所定範囲に収束するようにする手順を具備する請求項11に記載のX線CT画像再構成プログラム。
【請求項13】
前記各物質のX線吸収係数は予め設定された請求項9に記載のX線CT画像再構成プログラム。
【請求項14】
さらに、
前記物質kの係数akを
【数1】
を用いて演算する手順と、
該演算された係数ak及び前記発生X線連続エネルギースペクトルのエネルギーEiを用いて前記物質kのX線吸収係数μikを
【数2】
により演算する手順と
を具備する請求項9に記載のX線CT画像再構成プログラム。
【請求項15】
さらに、
前記物質kの係数akを
ak=αk・THn(n=1,2,…)+βk
但し、THnは前記しきい値、αk、βkは定数
により演算する手順と、
該演算された係数ak及び前記発生X線連続エネルギースペクトルのエネルギーEiを用いて前記各物質kのX線吸収係数μikを
【数2】
により演算する手順と
を具備する請求項7に記載のX線CT画像再構成プログラム。
【請求項16】
さらに、前記係数akにより前記物質kを推定する段階を具備する請求項14あるいは15に記載のX線CT画像再構成プログラム。
【図1】
【図2】
【図8】
【図10】
【図16】
【図20】
【図24】
【図25】
【図26】
【図31】
【図38】
【図40】
【図3】
【図4】
【図5】
【図6】
【図7】
【図9】
【図11】
【図12】
【図13】
【図14】
【図15】
【図17】
【図18】
【図19】
【図21】
【図22】
【図23】
【図27】
【図28】
【図29】
【図30】
【図32】
【図33】
【図34】
【図35】
【図36】
【図37】
【図39】
【図41】
【図42】
【図43】
【図2】
【図8】
【図10】
【図16】
【図20】
【図24】
【図25】
【図26】
【図31】
【図38】
【図40】
【図3】
【図4】
【図5】
【図6】
【図7】
【図9】
【図11】
【図12】
【図13】
【図14】
【図15】
【図17】
【図18】
【図19】
【図21】
【図22】
【図23】
【図27】
【図28】
【図29】
【図30】
【図32】
【図33】
【図34】
【図35】
【図36】
【図37】
【図39】
【図41】
【図42】
【図43】
【公開番号】特開2011−203160(P2011−203160A)
【公開日】平成23年10月13日(2011.10.13)
【国際特許分類】
【出願番号】特願2010−71819(P2010−71819)
【出願日】平成22年3月26日(2010.3.26)
【新規性喪失の例外の表示】特許法第30条第1項適用申請有り 平成21年9月28日 社団法人 計測自動制御学会計測部門発行の「第26回センシングフォーラム −センシング技術の新たな展開と融合− 資料」に発表
【出願人】(304021417)国立大学法人東京工業大学 (1,821)
【Fターム(参考)】
【公開日】平成23年10月13日(2011.10.13)
【国際特許分類】
【出願日】平成22年3月26日(2010.3.26)
【新規性喪失の例外の表示】特許法第30条第1項適用申請有り 平成21年9月28日 社団法人 計測自動制御学会計測部門発行の「第26回センシングフォーラム −センシング技術の新たな展開と融合− 資料」に発表
【出願人】(304021417)国立大学法人東京工業大学 (1,821)
【Fターム(参考)】
[ Back to top ]