【新規性喪失の例外の表示】特許法第30条第2項適用申請有り 掲載年月日:平成27年12月18日 https://www.jstage.jst.go.jp/article/sicetr/51/12/51_836/_article/−char/ja/
【解決手段】X線CT画像再構成法は、初期の第1CT画像を用意する第1の工程(ST1−4、1−5)と、第1CT画像に対して、エネルギー情報を考慮した順投影計算を行って、第1投影データを生成する第2の工程(ST1−6)と、第1投影データにFBP法による単一波長用画像再構成アルゴリズムを適用して第2CT画像を得る第3の工程(ST1−7)と、第2CT画像とオリジナルCT画像を比較し、その差分を算出する第4の工程(ST1−8)と、差分を第1CT画像μ´(x,y)に加算して得られるCT画像を、第1CT画像として設定する第5の工程(ST1−5)を含む。第2〜第5の工程を繰り返してメタルアーチファクトが低減されたCT画像を得る。
【背景技術】
【0002】
X線CT装置は、被写体の断層画像を非破壊で再構成する装置であり、医療、産業分野において広く使用されている。X線CT装置では、被写体を挟みX線源および検出器が対向配置され、これらが、被写体を中心として旋回可能となっている。X線CT画像の取得は、各旋回位置において被写体に対して、通常の波長分布を有するX線を照射し、検出器から得られるX線の透過強度を表す投影データを収集する第1のステップと、収集した投影データをコンピュータ上で処理して、被写体の断層画像であるX線CT画像を再構成する第2のステップからなる。
【0003】
第2のステップにおける投影データから被写体の断層画像を再構成するためのX線CT画像再構成方法としては、一般的に、解析的に断層画像を求めるフィルタ補正逆投影法(FBP法)に基づく方法、および、繰り返しを利用して徐々に断層画像を生成していく逐次近似画像再構成法と呼ばれる代数学的方法が用いられている。いずれの方法においても、CT画像再構成演算においては、物質のX線吸収係数のエネルギー依存性を考慮していない。すなわち、X線吸収係数μが一定の単一エネルギースペクトルのX線照射(単一波長照射)を前提とした単一波長用画像再構成アルゴリズムを用いている。
【0004】
X線CT画像撮影において、被写体の中に金属が含まれていると、メタルアーチファクトと呼ばれる放射状の虚像が発生する。メタルアーチファクトが発生すると、再構成画像であるX線CT画像における金属およびその周辺の形状が不鮮明となり、診断や検査を困難にしてしまう。
【0005】
メタルアーチファクトの主な発生要因は、実際に使用されている連続X線のエネルギースペクトルを考慮した画像再構成計算が行われていないことに起因している。実際に得られる投影データは連続エネルギースペクトルのX線照射によって得られたデータであるにも拘らず、X線吸収係数が一定の単一エネルギースペクトルのX線照射によって得られたデータであるとの前提で、画像再構成計算を行っている。この前提は、人体組織などのような非金属の被写体の場合には、X線吸収係数のエネルギー依存度が低いので、大きな問題にはならない。しかし、X線吸収係数の高い素材はエネルギー依存性が高いので、誤差を含む投影データが計算により生成され、このような投影データが原因となって、計算される再構成画像にメタルアーチファクトが発生する。
【0006】
特許文献1、2、3には、メタルアーチファクトを低減するX線CT画像再生方法が提案されている。特許文献1に開示の方法では、投影データのシノグラムから高X線吸収係数の物質領域、例えば金属領域を抽出して仮想投影によりその物質領域の投影データを改変して修正投影データを演算し、この修正投影データを逆投影してX線CT画像を再構成している。特許文献2に開示の方法では、投影条件を変更して複数回のCT撮影を行って複数の投影データを取得し、X線連続スペクトルの種類を限定して、X線吸収量を求めている。
【0007】
特許文献3に開示のX線CT画像再構成方法は、実際に測定されるオジリナル投影デー
タから再構成されたCT画像の画素値(CT値)に閾値処理を施すことによって物質領域を抽出し、これに基づき、連続エネルギースペクトルを考慮した投影データを計算によって求めている。そして、これと、単一エネルギースペクトルのX線を照射したときのX線投影データとの差から補正量を求め、この補正量を用いてオリジナル投影データを補正して得られる補正投影データを用いて画像再構成を行う処理を繰り返して、メタルアーチファクトを低減している。
【0008】
一方、非特許文献1には、逐次近似画像再構成によるアーチファクト低減手法が提案されている。その他にも、非特許文献2などにおいて逐次計算過程でアーチファクト低減を試みる手法は様々提案されている。
【発明の概要】
【発明が解決しようとする課題】
【0011】
メタルアーチファクトを低減するために、さまざまな手法による補正が考案されているが、未だ、この問題は十分に解決されていない。特に、逐次近似画像再構成法を用いたメタルアーチファクト低減は、現状では十分な成果がでていない。
【0012】
また、従来におけるメタルアーチファクト低減手法は、実際の測定で得られたオリジナルの投影データを利用している。一般に、X線CT装置では、再構成画像は保存されるものの、内部データである投影データは容量が大きいので保存されない。このため、メタルアーチファクト低減手法が確立されたとしても、投影データを利用する限り、そのような手法は、投影データが保存されていない過去のCT画像データに適用できない。
【0013】
本発明の課題は、測定されたオリジナルの投影データを必要とすることなく、当該投影データから再構成されたX線CT画像からメタルアーチファクトを低減可能なX線CT画像再構成方法およびX線CT画像再構成用プログラムを提案することにある。
【0014】
また、本発明の課題は、メタルアーチファクトを低減可能な逐次近似画像再構成法をベースにした新しいX線CT画像再構成方法およびX線CT画像再構成用プログラムを提案することにある。
【課題を解決するための手段】
【0015】
上記の課題を解決するために、本発明のX線CT画像再構成方法およびX線CT画像再構成用コンピュータプログラムは、
メタルアーチファクトが発生しているオリジナルCT画像を取得する第1の工程と、
各画素値を所定の値に設定した画像を第1CT画像として初期設定する第2の工程と、
前記第1CT画像に対して、連続エネルギースペクトルのX線照射における被写体を構
成する物質のX線吸収係数のエネルギー依存性を考慮した順投影計算を行って、第1投影データを生成する第3の工程と、
前記第1投影データに、単一エネルギースペクトルのX線照射を前提とする単一波長用画像再構成アルゴリズムを用いて逆投影計算を行って、第2CT画像を得る第4の工程と、
前記第2CT画像と前記オリジナルCT画像を比較し、その差分を算出する第5の工程と、
前記差分をスケーリングして前記第1CT画像に加算して得られるCT画像を、前記第1CT画像として再設定する第6の工程と
を有しており、
前記差分が所定値より小さくなるまで前記第3〜第6の工程を繰り返して、メタルアーチファクトが低減された前記オリジナルCT画像に対応する前記第1CT画像を得ることを特徴としている。
【0016】
本発明では、撮影で得られるオリジナルの投影データを使用せず、当該投影データから再構成されたオリジナルCT画像のみから、当該CT画像に発生しているメタルアーチファクトを低減することが可能である。メタルアーチファクト低減のために、被写体撮影で測定された膨大な容量の投影データを蓄積しておく必要がない。また、投影データが蓄積されていない過去のオリジナルCT画像に対してメタルアーチファクト低減処理を行うことができる。さらには、メタルアーチファクト低減のための計算処理によって、オリジナルCT画像の生成のために用いられた撮影で得られたオリジナルの投影データに対応する投影データも得ることができる。
【0017】
ここで、前記差分のスケーリングを適切に行うことにより、前記の第3〜第6の工程を繰り返す逐次計算における収束速度を速め発散の危険性を低減することができる。本発明者によれば、差分を5〜50分の1にスケーリングすれば良好な結果が得られることが確認された。
【0018】
次に、本発明のX線CT画像再構成方法およびX線CT画像再構成用コンピュータプログラムは、
金属を含む被写体に連続エネルギースペクトルのX線を照射して得られる当該被写体の投影データをオリジナル投影データとし、
前記オリジナル投影データに、単一エネルギースペクトルのX線照射を前提とする単一波長用画像再構成アルゴリズムを用いて逆投影計算を行って生成されるCT画像をオリジナルCT画像とすると、
メタルアーチファクトが発生している前記オリジナルCT画像の各画素値に基づき、当該オリジナルCT画像における金属領域を抽出し、当該金属領域を第1CT画像として初期設定する第1の工程と、
前記第1CT画像に対して逐次近似画像再構成法を適用して、当該第1CT画像を再設定する第2の工程と
を実行するようにしている。
【0019】
また、前記第2の工程では、前記第1CT画像に、前記連続エネルギースペクトルのX線照射における物質のX線吸収係数のエネルギー依存性を考慮した順投影計算を行って、第1投影データを生成し、前記第1投影データと前記オリジナル投影データを比較し、それらの差あるいは比を用いて、前記第1CT画像に対して逆投影計算を行って当該第1CT画像を更新する。
【0020】
そして、物質のX線吸収係数のエネルギー依存性を考慮した前記順投影計算においては、金属のX線吸収係数を用いて模擬した前記第1CT画像の画像濃度に対応するX線吸収
係数を使用し、前記第2の工程を1回あるいは複数回繰り返して、アーチファクトが低減された前記第1CT画像を得るようにしている。
【0021】
ここで、被写体を構成する物質のX線吸収係数のエネルギー依存性を考慮した前記順投影計算においては、前記第1CT画像の画像濃度をv、濃度諧調の数をn、画像濃度vに対応するX線吸収係数をμ
v、金属のX線吸収係数をμ
Metal、濃度分布の偏りに対応する係数をα(>1)とすると、次の式により規定されるX線吸収係数μ
vを使用して、前記第2の工程を1回あるいは複数回繰り返して、アーチファクトが低減された前記第1CT画像を得ることができる。
【数1】
【0022】
本発明者は、逐次近似画像再構成法における順投影計算において、上記のX線吸収係数分布の計算式を用いることにより、オリジナルX線CT画像に発生しているメタルアーチファクトを低減できることを確認した。
【0023】
一方、本発明では、前記オリジナルCT画像、前記第1CT画像、前記第2CT画像のうちの少なくとも一つの画像に対して二次元フーリエ変換を適用して、当該画像の低周波領域のパワースペクトルの総和SLP(Sum of Low frequency Power Spectrum)を求め、当該SLPの値を、前記画像におけるメタルアーチファクト発生量の定量的評価指標として用いている。
【0024】
数値シミュレーションにおける画像の評価にはRMSE(Root Mean Square Errors)が広く用いられているが、RMSEは基準となる画像が必要であるので、実際に撮影で得られるCT画像データに対しては利用できない。本発明者は、SLPとRMSEとの間に非常に高い相関があることに着目し、SLPを用いてX線CT画像を定量的評価指標として用いることを見出した。メタルアーチファクトの発生量を表す定量的指標としてのSLPを用いて、CT撮影における様々なパラメータを推定することが可能になり、また、上記のX線吸収係数分布の計算式に用いる最適パラメータを推定することが可能になる。
【発明を実施するための形態】
【0026】
[実施の形態1]
{X線CT装置の概要}
図1.1は、本発明の実施の形態1に係るX線CT装置の概略構成図である。X線CT装置1は、測定部2と、X線CT画像の画像再構成部3を含む制御部4と、入出力部5とを備えている。測定部2は、被写体6を挟み対向配置されるX線発生器7およびX線検出器8を備えている。
【0027】
撮影に当たっては、被写体6を中心としてX線発生器7およびX線検出器8を旋回させ、各旋回位置におけるX線照射により、X線検出器8からX線の透過強度を表す投影データが得られる。投影データは、コンピュータを中心に構成される制御部4の画像再構成部3において処理されて、被写体の断層画像であるX線CT画像が再構成される。再構成されたCT画像などの情報は、入出力部5の出力部、例えば表示部画面上に表示出力され、制御部4の内蔵メモリなどに記憶保持される。
【0028】
画像再構成部3では、X線検出器8から得られる投影データを、単一のエネルギースペクトルのX線を前提とした単一波長用画像再構成アルゴリズムを用いて逆投影計算を行って画像再構成を行う。この場合、入射X線強度をI
in、透過X線強度をI
out、被写体のある断面における物質のX線吸収係数の分布をμ(x,y)とすると、X線吸収係数とX線強度の関係は次の式(1.1)によって表される。ここで、X−Y座標系はx−y座標系を角度θだけ回転させた座標系である。
【0030】
入射X線強度I
inは既知であるので、投影データpは式(1.1)を変形して次の式(1.2)のように表される。
【数1.2】
【0031】
式(1.2)はμ(x,y)を、Y軸に沿って線積分した形になっており、この変形はラドン変換と呼ばれる。一方、逆投影とは、回転軸方向に向けて直前状に投影データを加えていく操作であり、その結果得られるCT画像g(x,y)は次の式(1.3)で表される。
【数1.3】
【0032】
なお、単純な逆投影では、再構成画像がボケてしまうので、一般的には、投影データに対してフィルタ処理を施してから逆投影を行うFBP法や、逐次近似法が用いられる。本例の画像再構成部3では、FBP法による単一波長用画像再構成アルゴリズムを用いて逆投影演算を行ってCT画像を再構成している。
【0033】
ここで、実際に使用されているX線源は連続エネルギースペクトルのX線である。このような連続X線では、式(1.1)〜(1.3)で表される関係が成立しない。物質のX線吸収係数μは入射X線のエネルギーに依存するので、式(1.1)は厳密には次の式(1.4)に示すエネルギー積分系で表される。
【0035】
このとき、投影データpは、次の式(1.5)で表される。
【数1.5】
【0036】
しかし、エネルギー別のX線強度を取得するのは困難であるので、多くのX線CT装置と同様に、本例においても、画像再構成部3では、式(1.2)および式(1.3)を用いた画像の再構成計算を行っている。多くの場合、この近似計算に伴った問題は生じないが、被写体に金属のようなX線吸収係数のエネルギー依存性が高い物質が含まれていると、ビームハードニング(X線スペクトルが高エネルギー側に集中する現象)が生じ、投影データの物質厚みに対する線形性が損なわれる。投影データの非線形性が計算に矛盾を発生させ、この矛盾がメタルアーチファクトの原因となる。
【0037】
そこで、本例の画像再構成部3には、メタルアーチファクト低減部10が含まれている。メタルアーチファクト低減部10は、以下に述べるように、エネルギー情報を考慮して、測定された投影データに基づき計算されたCT画像に含まれるメタルアーチファクトの低減処理を行い、メタルアーチファクトが除去あるいは低減されたCT画像を生成する。
【0038】
{メタルアーチファクト低減方法}
図1.2は、主としてX線CT装置1のメタルアーチファクト低減部10において行われるメタルアーチファクト低減の手順を示す概略ブロック図である。まず、X線CT装置1は、被写体撮影を行って(ST1−1)、X線検出器8から得られる測定投影データであるオリジナル投影データp
originalを、一般的に用いられているフィルタ補正逆投影法(FBP法)、逐次近似法等の単一波長用画像再構成アルゴリズムを用いて逆投影計算を行って、CT画像を再構成する(ST1−2、ST1−3)。例えば、FBP法を用いて逆投影計算を行う。被写体に金属等が含まれている場合には再構成されたCT画像にはメタルアーチファクトが発生する。被写体断面におけるX線吸収係数の分布をμ(x,y)、X−Y座標系をx−y座標系を角度θ(θ:X線投影角度)だけ回転させた座標系とし、メタルアーチファクトが発生しているCT画像をオリジナル再構成画像μ
original(x,y)とする。
【0039】
画像再構成部3のメタルアーチファクト低減部10において、オリジナルCT画像μ
original(x,y)からのメタルアーチファクト低減処理は、次の手順によって行われる。
【0040】
手順1:初期の画像を用意し、これを第1CT画像μ´(x,y)として初期設定する(ST1−4、ST1−5)。初期の第1CT画像μ´としては、例えば、オリジナルCT画像の画素値の平均値を一様に分布させた画像を用いる。
【0041】
手順2:初期の第1CT画像μ´(x,y)に対して、被写体を中心とする所定のX線投影角度θの方向から、エネルギー情報を考慮した順投影計算を行い、式(1.6)で表される第1投影データp´(X,θ)を得る(ST1−6)。
【0043】
手順3:手順2の結果得られた第1投影データp´(X,θ)にFBP法を適用し、式(
1.7)で表される第2CT画像μ”(x,y)を得る(ST1−7)。この第2CT画像の計算は、画像再構成部3においてX線検出器8から得られるオリジナル投影データp
originalからオリジナルCT画像μ
original(x,y)を計算する場合と同様に、エネルギー情報を考慮せずに行う。本例ではFBP法による単一波長用画像再構成アルゴリズムを用いて逆投影計算を行って、第2CT画像を再構成する。
【0045】
手順4:第2CT画像μ”(x,y)とオリジナルCT画像μ
original(x,y)を比較し、式(1.8)に示すように、その差分dμを算出する(ST1−8)。
【数1.8】
【0046】
手順5:差分をフィードバック係数wでスケーリングし、これを、初期の第1CT画像μ´(x,y)に加算し、得られたCT画像を、式(1.9)に示すように、初期の第1CT画像の代わりに、第1CT画像として再設定する(ST1−5)。
【0048】
手順6:手順2から手順5をすべての投影角度θについて、フィードバック値(差分dμ)が十分に小さくなるまで、繰り返す(
図1.2のST1−5〜ST1−8)。式(1.10)、(1.11)に、k回(kは正の整数)の繰り返しによって得られる第1投影データp
kおよび、当該第1投影データp
kを用いて更新された第1CT画像μ
k+1をそれぞれ表す。これにより、実際のX線吸収分布に対応するアーチファクトが十分に低減された第1CT画像μ
k+1が得られる。
【0050】
このように、実施の形態1に係るアーチファクト低減方法では、オリジナル投影データp
originalを必要とせず、オリジナルCT画像μ
originalのみがあればよい。本方法では、アーチファクトが発生しているオリジナルCT画像μ
originalが算出される、元となるオリジナル投影データp
originalをシミュレーションによって、第1投影データp
kとして計算により算出している。すなわち、エネルギー情報を考慮せずに画像再構成を行って得られるメタルアーチファクトが発生しているCT画像が、どのような投影データ(第1投影データ)から得られるのかも計算によって算出できる。
【0051】
{エネルギー情報の考慮:X線吸収係数の仮定}
エネルギー情報を考慮した順投影計算(上記の手順2:ST1−6)を行うためには、CT画像における各物質の画素値(CT値)、したがって、各物質のエネルギーEに対応するX線吸収係数μ(E)の情報が必要になる。本例では、5000種類の物質について−1000から4000HUの範囲のX線吸収係数分布を用いた。この分布は、表1.1に示す6種類の物質のX線吸収係数から線形補間を行って算出したものである。なお、CT画像の画素値が−1000を下回る場合あるいは4000を上回る場合には、画素値をそれぞれ−1000、4000に修正して用いた。
【0053】
{シミュレーション実験}
実施の形態1によるメタルアーチファクト低減方法の効果を確認するためにシミュレーション実験を行った。
【0054】
(実験条件)
Shepp−Loganファントム(以下、単に「ファントム」と呼ぶ場合もある。)に金属を埋め込んだモデルを用意して、フィードバック係数α=10でのシミュレーション実験を行った。以下、すべての画像において、画像サイズ=512×512、FOV=300mm、WW=1500、WL=300である。
【0055】
本実験では、実際の状況を想定して、FBP結果のCT画像を16bit(−32768〜32767)の範囲に丸めている。また、実際のCT値と対応させるために、空気を−1000、脂肪を−100、水(軟組織)を0、硬組織を150、骨を1000、金属(鉄)を4000としており、これらのエネルギーに対する質量吸収係数分布はNISTのデータベースより引用した。計算過程で必要な、これらの中間を含む合計5000のX線吸収係数分布は、各物性値を線形補間により算出している。なお、WW=1500、WL=300で観察できるCT値の範囲は−450〜1050であり、結果のCT画像上では骨と金属は色の区別がほとんどつけられない。
【0056】
メタルアーチファクト低減の評価指標としてRMSE(2乗平均平方根誤差)を用いた。以下において、RMSE値は、ファントムと、実施の形態1による方法による再構成画像(FBP画像あるいはFBP値)との間の誤差を表している。RMSE
m値は、メタルアーチファクトの発生している元のオリジナルCT画像と、逐次計算の過程で得られるCT画像(第2CT画像)との間の誤差を表している。なお、メタルアーチファクトの発生量の定量的指標として、RMSEの代わりに、後述のパワースペクトルの総和SLP(Sum of Low frequency Power Spectrum)を用いることが可能である。
【0057】
(実験その1)
図1.3(a)は、金属(鉄)が埋め込まれていない基本のShepp−Loganファントムを示し、
図1.3(b)はそのFBP法による結果を示す。
図1.3(c)は6個の金属が埋め込まれたファントムを示し、
図1.3(d)にはそのFBP画像(RMSE:666.47)を示す。また、
図1.4は、逐次計算の過程で得られた再構成画像(繰り返し回数:1、10、30、100、500、1000)を示し、
図1.5は逐次計算に伴うRMSEの変化を示すグラフである。
【0058】
これらの図からわかるように、元のオリジナルCT画像との差が激しい部分から優先的に逐次過程のCT画像に反映されている様子が伺える。ある程度、画像が生成されてくると、元のCT画像と等価な画像を生成するために、細かい部分の調整が少しずつ行われていく。最終的に、RMSEは22.775まで、RMSE
mは1.4920まで減少しており、僅かにメタルアーチファクトの筋が残っているが、定性的にも有用な画像であると判断できる。
【0059】
(フィードバック係数w)
ここで、フィードバック係数wの値を小さくすれば、収束速度が速くなるが、画像が収束しない可能性が高くなる。また、収束しているように見えても微小なCT値の反映が困難であるため、ノイズが残りやすい傾向にある。逆にwの値が大きいと収束は遅くなるが、発散の危険性が大幅に減る。ただし、あまり大きすぎると、ローカルミニマムに陥り、抜け出せなくなる可能性があると考えられる。本シミュレーションにおいては、少なくとも、
5 ≦ w ≦ 50
の範囲内であれば、最終的に同等の結果が得られることが分かった。なお、計算時間は1回の逐次計算で20秒、1000回の計算で6時間程度であった。
【0060】
(実験その2)
より複雑なファントムを用意して効果の検証を行った。上記の(実験その1)で使用した6個の金属が埋め込まれたファントムに、様々な形状の金属を埋め込んだモデルを用いた。
図1.6には、ファントム、そのFBP結果、逐次計算の過程で得られた再構成画像(繰り返し回数:30、100、1000)をそれぞれ示す。結果はいずれも、WW=1500、WL=300におけるものである。
【0061】
FBP画像で激しく発生していたメタルアーチファクトは、1000回の繰り返しで大幅に低減されていることが確認できる。複雑な金属の形状も、ほぼ忠実に再構成されている。凹凸の部分アーチファクトが強く残っている領域が何か所か見られる。RMSE
mの値が小さくなるにつれて、画像の更新量も減るため、メタルアーチファクト低減の速度は分数関数的に減っていく。更新回数を大幅に増やせば、ここで残ったアーチファクトも低減される可能性がある。
【0062】
{その他の実施の形態}
図1.1、1.2に示す例は、X線CT装置1に本発明によるアーチファクト低減方法を適用したものである。本発明は、既に得られているオリジナルCT画像を取り込み、このオリジナルCT画像からメタルアーチファクトが低減されたCT画像を再構成するCT画像再構成用の処理装置、処理方法あるいは処理プログラムとして用いることもできる。この場合には、X線発生器、X線検出器が備わっていない装置構成とされる。
【0063】
[実施の形態2]
{X線CT装置の概要}
図2.1は本発明の実施の形態2に係るX線CT装置の概略構成図である。X線CT装置11の基本構成は、実施の形態1のX線CT装置と同様であり、測定部12と、X線CT画像の画像再構成部13を含む制御部14と、入出力部15とを備えている。測定部12は、被写体16を挟み対向配置されるX線発生器17およびX線検出器18を備えている。
【0064】
画像再構成部13では、X線検出器18から得られる投影データを、単一のエネルギースペクトルのX線を前提とした単一波長用画像再構成アルゴリズムを用いて逆投影計算を行って画像再構成を行う。すなわち、実際に使用されているX線源は連続エネルギースペクトルのX線であるが、エネルギー別のX線強度を取得するのは困難であるので、多くのX線CT装置と同様に、本例においても、画像再構成部13では、単一のエネルギースペクトルのX線を前提とした単一波長用画像再構成アルゴリズムを用いて逆投影計算を行って画像再構成を行う(前述の式(1.2)および式(1.3)参照)。
【0065】
したがって、X線CT装置11においても、被写体に金属のようなX線吸収係数のエネルギー依存性が高い物質が含まれていると、再構成画像にメタルアーチファクトが発生する場合がある。
【0066】
そこで、画像再構成部13には、メタルアーチファクト低減部20が含まれている。メタルアーチファクト低減部20は、以下に述べるように、エネルギー情報を考慮して、逐次近似画像再構成手法を用いて、測定された投影データに基づき計算されたCT画像に含まれるメタルアーチファクトの低減処理を行い、メタルアーチファクトが低減されたCT画像を再構成している。また、本例では逐次近似画像再構成手法として、ART(Algebraic Reconstruction Technique)法を採用している。
【0067】
{逐次近似画像再構成手法を利用したメタルアーチファクト低減方法}
図2.2は、画像再構成部13における逐次近似画像再構成手法を利用したメタルアーチファクト低減方法の流れ図である。金属を含む被写体16に連続エネルギースペクトル
のX線を照射して得られる当該被写体の投影データをオリジナル投影データと呼び、このオリジナル投影データに、単一エネルギースペクトルのX線照射を前提とするFBP法による単一波長用画像再構成アルゴリズムを用いて逆投影計算を行って生成されるCT画像をオリジナルCT画像と呼ぶものとすると、本例の方法は次の2つの工程からなる。
【0068】
<第1の工程>
FBP法適用後の閾値処理、もしくは、先に述べた特許文献2で提案している手法によって、計測された投影データからFBP法による単一波長用画像再構成アルゴリズムを適用して再構成したオリジナルCT画像から金属領域を抽出し、これを、逐次計算における初期のCT画像である第1CT画像として設定する(ST2−11、ST2−12)。
<第2の工程>
エネルギー依存性を考慮した順投影計算を順問題とした逐次近似画像再構成法であるART法を適用して、第1CT画像を逐次更新して、メタルアーチファクトが低減された第1CT画像を演算により求める(ST2−12、ST2−13、ST2−14)。
【0069】
具体的に説明すると、第1の工程に先立って、X線発生器17から連続エネルギースペクトルを有するX線を被写体16に照射し(ST2−1)、画像再構成部13において、X線検出器18によって計測されたオリジナル投影データを収集し(ST2−2)、収集したオリジナル投影データから、FBP法による単一波長用画像再構成アルゴリズムを適用してオリジナルCT画像を再構成する(ST2−3)。第1の工程では、アーチファクト低減部20において、再構成したオリジナルCT画像から、その各画素値に基づき、金属領域を抽出する(ST2−11)。抽出した金属領域を、メタルアーチファクト低減処理を行うための繰り返し計算のための初期の第1CT画像として用いる(ST2−12)。
【0070】
FBP法による画像再構成では、メタルアーチファクトによりコントラスト分解能が低下し、非金属領域の観察が困難になる。金属領域は顕著な値を示すため、その抽出は画素値の閾値処理によって精度よく行えることが多い。金属領域の抽出が困難な場合には、特許文献2に記載の手法が有用である。すなわち、オリジナルCT画像の画素値を基準に画素数のヒストグラムを作成し、これに基づき閾値を設定して、金属領域を抽出する。
【0071】
次に、第2の工程では、アーチファクト低減部20において、エネルギー依存性を考慮した順投影計算を順問題とした逐次近似画像再構成法であるART法を適用して、メタルアーチファクトが低減された推定CT画像を演算により求める(ST2−12〜ST2−14)。順投影計算の際にX線吸収係数のエネルギー依存性を考慮することで、逐次計算における矛盾を解消している。単純な差や積のフィードバックでは、金属領域と非金属領域の差を明確に識別することは困難であるが、あらかじめ抽出した金属領域を初期の第1CT画像とすることで、高速で正確な収束が期待される。なお、本手法では、金属領域の抽出が正確に行えることを前提とし、逐次近似計算において金属領域のフィードバックは行わないものとした。
【0072】
さらに具体的に説明すると、逐次近似画像再構成法における逆投影(フィードバック)の種類には加算型と乗算型が存在し、加算型では投影データの差を、乗算型では投影データの比を逆投影する。本例では、逐次近似画像再構成法であるART法として加算型ART法を採用している。加算型ART法は次の手順からなる。
手順1:初期の第1CT画像μ´(x,y)を用意する(ST2−12)。
手順2:初期の第1CT画像μ´(x,y)に対して、エネルギー依存性を考慮して、投影角度θの方向からの順投影処理を行う(順投影計算)。この結果、X線の第1投影データp´(X,θ)を得る(ST2−13)。
手順3:手順2の結果得られた第1投影データp´(X,θ)と、測定部12による実際
の計測で得られているオリジナル投影データp(X,θ)を比較し、その差あるいは比をとり、これを初期の第1CT画像μ´(x,y)にフィードバック(逆投影計算)する(ST2−14)。
手順4:手順2および手順3をすべての投影角度θに対して行う。
手順5:手順2〜4を1回、あるいは複数回繰り返す。一般的には、フィードバック量(差あるいは比)が十分に少なくなるまで繰り返す。これにより、実際のX線吸収分布を表すCT画像μ(x,y)を得る。
【0073】
X線吸収係数のエネルギー依存性を考慮したART法の繰り返しの式は、k回目の第1CT画像をμ(x、y、E)
(k)、k+1回目の第1CT画像をμ(x、y、E)
(k+1)とすると、式(2.1)で表される。
【0075】
以下において、エネルギー情報を考慮した本例のART法のことを、接頭辞「E」をつけて、「E−ART法」と呼ぶ。
【0076】
{エネルギー情報の考慮:X線吸収係数の仮定}
エネルギー情報を考慮した順投影計算(上記の手順2:ST2−12、ST2−13)を行うためには、各物質のエネルギーEに対応するX線吸収係数μ(E)の情報が必要になる。エネルギーに対する吸収係数分布はNISTのデータベースに基づき算出した。実際の元素のX線吸収係数をグラフ化したものを
図2.3に示す。
【0077】
ここで、逐次近似計算においては、逆投影画像の画像濃度が吸収係数を表していると捉えることができるので、逆投影画像の画像濃度vに対応するX線吸収係数μv(E)があればよい。そこで、実施の形態2では、金属のX線吸収係数μ
Metal(E)のみを用いて
図2.3のグラフを模擬したX線吸収係数μv(E)を得るために、次の式(2.2)を採用している。
【0079】
ここで、αは濃度分布の偏りに対応する係数、nは濃度階調の数を表す。αは1より大きい値であり、αが大きいほど非金属付近のX線吸収係数が密に割り当てられ、αが1に近いほど空気から金属にかけてのX線吸収係数が等間隔に割り当てられる。このように、割り当てられるX線吸収係数の密度を変化させることで、断層面に存在する物質の分布状況に応じた繰り返し計算が期待される。
【0080】
例えば、金属以外のX線高吸収物質が存在しない被写体の場合は、非金属付近のX線吸収係数を密に割り当てた方が、繰り返し計算の過程で実際の物質に近いX線吸収係数に収束する可能性が高くなる。
【0081】
本例では、離散スペクトルを持った離散X線に対する吸収係数を考え、μ
Metalとして、鉄のX線吸収係数(表2.1参照)を用いた。離散スペクトルは、スペクトルアナライザ(AMPTEX Inc.: XR-100T-CdTe)を用いて、マイクロフォーカスX線源(浜松ホトニクス(株):L9181−02)に管電圧130kVを加えたときのX線スペクトルから取得したものである。
【0083】
n=256として式(2.2)によりX線吸収係数μv(0 ≦ v ≦ 255)を
算出し、μ
16i(0 ≦ i ≦ 15)およびμ
255(=μ
Metal)をグラフ化した結果を
図2.4に示す。
【0084】
図2.4(a)はα=1.10、
図2.4(b)はα=1.05、
図2.4(c)はα=1.03となっており、X線吸収係数の分布はαの値が1.0から遠いほど偏りが大きく、αの値が1.0に近いほど等間隔になることが観察できる。また、
図2.3と比較することで、本手法により仮定したX線吸収係数は実際のX線吸収係数を模擬できると言える。
【0085】
図2.5には、n=256、α=1.03におけるX線吸収係数分布を利用したシミュレーション結果を示してある。
図2.5(a)の数値ファントムは、樹脂の中に鉄が埋め込まれた実験用サンプルを想定しており、樹脂部分のX線吸収係数は、画像濃度v=80のときのX線吸収係数μ
80(E)の値を設定している(
図2.4(c))。この数値ファントムに周囲512方向からの順投影を行い、FBP法による画像再構成を行った。入射X線として単色X線を設定したシミュレーション結果を
図2.5(b)に、入射X線として表2.1の離散スペクトルを設定したシミュレーション結果を
図2.5(c)に示す。
図2.5(d)にはシミュレーション条件を纏めて示してある。単色X線では計算に矛盾が生じないため、メタルアーチファクトが一切発生していない。一方で、離散X線を使用した場合はメタルアーチファクトが発生し、樹脂領域および金属周辺部が不鮮明になっていることが確認できる。
【0086】
{二次元フーリエ変換による画像の評価}
数値シミュレーションにおけるCT画像の評価には、RMSE(Root Mean Square Errors)が広く用いられているが、RMSEは基準となるCT画像が必要であるので、実際に撮影して得たオリジナルCT画像データに対しては利用できない。そこで、本例では、二次元フーリエ変換を用いたアーチファクト発生量の評価手法を提案する。アーチファクトの発生量を表す定量的指標があれば、CT撮影におけるさまざまなパラメータを推定することが可能になる。
【0087】
まず、得られた再構成画像に対して二次元フーリエ変換を適用し、低周波領域のパワースペクトルの総和(Sum of Low frequency Power Spectrum:SLP)を求める。画像をf(x,y)、画像を二次元フーリエ変換した結果をF(u,v)、F(u,v)の低周波領域のパワースペクトルの総和をSLP(F)とすると、それぞれの関係式は次の式(2.3)、(2.4)で表される。ここで、Cは低周波領域を表している。本例では、二次元フーリエ変換の結果に対して、画像の中心から画像一辺の1/20の半径を持つ円形領域を低周波領域と定義している。
【0089】
図2.6(a)は、
図2.4(b)の中央部のラインプロファイル(A−A´)、
図2.6(b)は
図2.4(c)の中央部のラインプロファイル(B−B´)である。
図2.6(b)を見ると、金属の両端で画素値のオーバーシュートが生じており、これによって非金属領域のコントラストが低下し、樹脂領域と空気領域の境界が不鮮明になっていることが分かる。この2つのラインプロファイルをフーリエ変換した結果を
図2.7に示す。
【0090】
この結果から、金属による画素値のオーバーシュートが少ない方が、すなわち、メタルアーチファクトの発生が少ない方が、低周波領域のパワースペクトルが増加する傾向にあることが分かる。これは、金属による信号のオーバーシュートによって信号全体が底上げされ、低周波成分による信号の再現が困難になるためであると考えられる。
【0091】
以上のことから、画像全体の低周波領域のパワースペクトルの総和であるSLP(F)は、メタルアーチファクトが少ないほど増加することが予想される。このことを確認するまえに、メタルアーチファクト発生量を調整した画像を幾つか生成し、SLPとRMSEとの相関分析を行った。
【0092】
メタルアーチファクト発生量の調整は、
図2.4(a)の金属領域の物性を変化させることで行っており、
図2.8(a)ではマグネシウム、
図2.8(b)ではアルミニウム、
図2.8(c)ではカルシウム、
図2.8(d)ではチタンをそれぞれ設定している。定性的にも、金属の原子番号が高くなるにつれて、メタルアーチファクト発生量が増加していることが確認できる。これらの画像に対するSLPと、真値である
図2.4(a)のファントムとのRMSEの関係をプロットした結果を、
図2.9に示す。
【0093】
決定係数はR
2=0.9328と非常に高い値を示しており、これは、SLPとRMSEとの間には、非常に強い相関がある、すなわち、SLPはアーチファクト発生量の評価に利用可能であることを意味している。
【0094】
{検証実験その1}
X線CT装置(コムスキャンテクノ(株):ScanXmate-A130SS940)を用いて、被写体撮影を行い、得られたデータに対して検証を行った。X線CT装置(
図2.1参照)の各部は次の通りである。
検出器:東芝電子管デバイス(株)、9インチI.I.システム、E5761SCA1−2NX線発生器:浜松ホトニクス(株)、反射型密閉管、L9181−02
最大管電圧130kV、最大管電流300μA、焦点寸法20μm
ただし、本実験では、小焦点モードを用いており、このときの最大管電圧は130kV、最大管電流は200μA、焦点寸法は8μmである。
【0095】
図2.10(a)に実験で用いた試験片を示す。この試験片は、アクリル製の樹脂とビーズ、および鉄製の釘で構成されており、CT撮影およびFBP法を適用すると、
図2.10(b)の画像が得られる。釘の部分からメタルアーチファクトが発生することで、周辺部の断面形状が不鮮明になっており、画像全体のコントラストも悪くなっていることがわかる。なお、
図2.10(c)には投影条件を纏めて示してある。
【0096】
図2.11は、
図2.10のデータに対してE−ART法を適用した際の、αに対するSLPの変動を示したグラフである。SLPの値は、α=1.036において、最大値8546.8をとっており、このとき最もメタルアーチファクトが低減されていることが予想される。
【0097】
図2.12に、各αにおけるE−ART法の画像再構成の結果を示す。α=1.010からα=1.036にかけて、メタルアーチファクトの大幅な低減が、α=1.036からα=1.100にかけて、メタルアーチファクトの増加が確認できる。このことから、SLPによるアーチファクトの定量的評価は、定性的な評価と一致しているといえる。
【0098】
{検証実験その2}
図2.12に示した結果は、いずれも、E−ART法を1回だけ繰り返したものとなっている。この中で最も良好な結果を示した、α=1.036におけるE−ART法を5回繰り返した結果を
図2.13(a)に示す。
図2.12(c)と同様の結果ではあるが、樹脂領域の輪郭が鮮明になり、SLPも8655.8に向上している。
【0099】
図2.13(b)は、
図2.10のファントムから釘を取り除いて撮影、再構成した結
果を示す。このファントムは金属が含まれていないので、メタルアーチファクトは一切発生していない。その代りに、リングアーチファクトが強調される結果となっているが、すべてのビーズの断面が正確に観察可能となっており、良好な再構成画像といえる。
図2.13(b)のSLPは8676.4となっており、この数値からもSLPによる評価の妥当性を確認することができる。
【0100】
以上説明したように、本実施の形態では、メタルアーチファクトを低減するために、順投影計算におけるX線吸収係数分布の計算式を構築し、エネルギー情報を考慮した新しい逐次近似画像再構成手法を提案した。また、X線吸収係数分布の計算式に用いる最適パラメータを推定するために、二次元フーリエ変換のパワースペクトルを用いた、アーチファクト発生の強弱を表す定量化指標を提案した。
【0101】
樹脂と金属で構成された実験用ファントムを用意し、CTで撮影したデータに対して本提案手法を適用することで、定量化指標の妥当性、および、手法の有効性が示された。
【0102】
実験結果から、実施の形態2に係るエネルギー情報を考慮した逐次近似法による画像再構成は、メタルアーチファクトの低減に効果的であることが確認できた。