JP3707347B2 - X線ct装置の画像処理方法及びx線ct装置並びにx線ct撮影用記録媒体 - Google Patents
X線ct装置の画像処理方法及びx線ct装置並びにx線ct撮影用記録媒体 Download PDFInfo
- Publication number
- JP3707347B2 JP3707347B2 JP2000106271A JP2000106271A JP3707347B2 JP 3707347 B2 JP3707347 B2 JP 3707347B2 JP 2000106271 A JP2000106271 A JP 2000106271A JP 2000106271 A JP2000106271 A JP 2000106271A JP 3707347 B2 JP3707347 B2 JP 3707347B2
- Authority
- JP
- Japan
- Prior art keywords
- image
- projection data
- estimated
- ray
- updated
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Lifetime
Links
- 238000003672 processing method Methods 0.000 title claims description 19
- 238000013170 computed tomography imaging Methods 0.000 title claims description 16
- 238000000034 method Methods 0.000 claims description 176
- 239000006096 absorbing agent Substances 0.000 claims description 59
- 238000009795 derivation Methods 0.000 claims description 23
- 238000012937 correction Methods 0.000 claims description 20
- 230000000052 comparative effect Effects 0.000 claims description 12
- 238000001914 filtration Methods 0.000 claims description 11
- 238000005259 measurement Methods 0.000 claims description 8
- 230000007423 decrease Effects 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 description 21
- 230000006870 function Effects 0.000 description 21
- 230000000694 effects Effects 0.000 description 15
- 230000005540 biological transmission Effects 0.000 description 12
- 238000010586 diagram Methods 0.000 description 12
- 238000012545 processing Methods 0.000 description 11
- 238000001514 detection method Methods 0.000 description 9
- 238000010521 absorption reaction Methods 0.000 description 8
- 238000003384 imaging method Methods 0.000 description 6
- 238000013480 data collection Methods 0.000 description 5
- 238000004422 calculation algorithm Methods 0.000 description 3
- 230000002093 peripheral effect Effects 0.000 description 3
- 239000012466 permeate Substances 0.000 description 3
- 230000003252 repetitive effect Effects 0.000 description 3
- 230000002745 absorbent Effects 0.000 description 2
- 239000002250 absorbent Substances 0.000 description 2
- 210000000746 body region Anatomy 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 239000002184 metal Substances 0.000 description 2
- 239000000047 product Substances 0.000 description 2
- 230000012447 hatching Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y10—TECHNICAL SUBJECTS COVERED BY FORMER USPC
- Y10S—TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y10S378/00—X-ray or gamma ray systems or devices
- Y10S378/901—Computer tomography program or processor
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Description
【発明の属する技術分野】
この発明は、X線CT装置の画像処理方法及びX線CT装置並びにX線CT撮影用記録媒体に係り、特に、金属等のX線高吸収体によるX線の吸収または散乱によって断層画像上に発生する偽像を低減させる技術に関する。
【0002】
【従来の技術】
一般にX線CT装置は、図8に示すように、X線ビームBを発生させるX線管51と、X線を検知するX線検知器52とを備えている。そしてX線管51とX線検知器52とは、被検体Pを挟んで備えられており、被検体Pの体軸Z(図8の紙面に対して垂直に延びる方向)の周りに回転しながら、X線管51は被検体PにX線ビームBを曝射して撮影を行う。一般にX線検知器52は、複数個(この場合n個)の検出素子CHi (i=1,2,.......,n-1,n nは自然数) から構成されており、それらの検出素子CHi は微小角毎に、扇状に配列されて設置されている(図8中の▲1▼)。
【0003】
従来の画像処理方法の具体例では、図8に示すように、被検体Pの周りからX線を曝射して、投影平面F上にX線透過データf(P)を得る(図8中の▲2▼)。なお、本明細書では、これを「被検体PをX線ビームBの曝射方向に投影する」(以下、適宜「被検体Pを投影する」と略記する)ともいい、このとき得られる投影平面F上のX線透過データf(P)を実測投影データf(P)ともいう。さらに投影平面F上の前記実測投影データf(P)にフィルタリング及び逆投影等の再構成処理を施して、断層平面G上に実測断層画像G[f(P)] を得る(図8中の▲3▼)。本明細書では、コンボリューション(Convolution) 関数を用いて畳み込み積分(重畳積分)を行なうことを、フィルタリングと定義付ける。また、投影データにフィルタリング及び逆投影の再構成処理を施す方法は、一般にFBP(Filtered Back Projection)法として知られている。なお、FBP法は代表的な再構成処理方法の1つなので、その説明については省略する。なお、ある断層画像を数学的アルゴリズムを用いて投影データを導出することを、「断層画像をX線の曝射方向に順投影する」(以下、適宜「断層画像を順投影する」と略記する)という。
【0004】
なお投影平面Fの横軸方向は検出素子CHi の並びを示し、縦軸方向はX線ビームBの曝射方向の角度θの並びを示す。また、投影平面F及び断層平面G内のハッチング領域は模式図である。また、投影平面F及び断層平面Gは領域、即ち各々の座標の集合を表して、各座標には数値の重みをそれぞれ持たないものとして、以下の説明を行う。
【0005】
なお、ある物体がPで表されるとき、その物体Pの投影はf(P)で表して、断層画像がαで表されるとき、断層画像αの順投影はF(α)で表して、ある投影データがβで表されるとき、その投影データβのフィルタリングと逆投影とによって得られた断層画像、即ち、FBP法によって再構成された断層画像はB[H(β)] またはG(β) で表すとする。このとき、逆投影の演算子はBで表し、フィルタリングの演算子はHで表し、逆投影とフィルタリングとによる演算子はGで表すとする。また、被検体の投影と断層画像の順投影とを区別するためにそれぞれの投影と順投影とに関する演算子は、fとFとで表すとする。また、測定誤差があるので、投影データβが実測のとき、その実測投影データβのFBP法によって得られた断層画像を順投影しても、元の実測投影データβと全く同じデータには戻らない。従って、本明細書ではβ=F[G(β)]は成立しないものとして、以下の説明を行う。なお、投影データ(実測投影データf(P)も含む)及び断層画像(実測断層画像G[f(P)] も含む)は、投影平面F及び断層平面Gと違い、数値の重みを持つもの、即ち本明細書では画素値であるとして、以下の説明を行う。
【0006】
【発明が解決しようとする課題】
しかしながら、従来の画像処理方法の場合には、次のような問題がある。
即ち、被検体を撮影して断層画像を再構成すると、偽像が発生する。特に、金属等に代表される高吸収体を含む被検体を撮影する場合、高吸収体の吸収または散乱によって、高吸収体及び高吸収体周辺部に発生する偽像は顕著なものになる。以下、本明細書では上記のような偽像をアーティファクトという。その中でも特に、高吸収体を中心として放射状の縞模様が発生するストリークアーティファクトと、複数個の高吸収体に挟まれた部分に生じる影状のシェーディングアーティファクトとが知られている。それらを含んだアーティファクトを低減させるために、様々な手法が提案されている。代表的な低減手法としては、下記の2つの手法、即ちIRR(Iterative Reconstruction and Reprojection) 法と、ART(Algebraic Reconstruction Technique)/EM(Expectation and Maximization)法とが知られている。
【0007】
先ず、IRR(Iterative Reconstruction and Reprojection) 法について、図9を参照して、説明をする。上記のIRR法では、従来の画像処理方法で得られた断層平面G上の実測断層画像G[f(P)] を再度X線の曝射方向に数学的アルゴリズムによって順投影して、投影平面F上に修正用の順投影データF(G[f(P)])を得る。そしてオペレータは断層平面G上の実測断層画像G[f(P)] を参照しながら、高吸収体に相当する部分を高吸収体領域Mとして断層平面G上に設定する(図9中の▲1▼)。なお、高吸収体領域Mは閉じられた領域である。そして、高吸収体領域を透過した部分Lが投影平面F上に作成される。高吸収体領域を透過した部分Lに関して実測投影データf(P)の画素値を、修正用の順投影データF(G[f(P)])の画素値或いは修正用の順投影データF(G[f(P)])に応じた画素値に置換して、実測投影データf(P)を修正して修正投影データF(P1) を求める(修正用の順投影データF(G[f(P)])と修正投影データF(P1) とは別異であることに注意されたい)(図9中の▲2▼)。さらにその修正投影データF(P1) の再構成を再度行って、断層平面G上に修正断層画像G[F(P1)]を得る(図9中の▲3▼)。
【0008】
なお、高吸収体領域を透過した部分Lに関して実測投影データf(P)の画素値を、修正用の順投影データF(G[f(P)])の画素値に置換するとは、詳述すると修正投影データF(P1) を下記の(1)式及び(2)式のように求めることである。
*高吸収体領域を透過しなかった部分では
F(P1) =f(P) ・・・(1)
*高吸収体領域を透過した部分Lでは
F(P1) =F(G[f(P)]) ・・・(2)
なお、高吸収体領域を透過した部分Lに関して実測投影データf(P)の画素値を、修正用の順投影データF(G[f(P)])に応じた画素値に置換する一例として、下記の(3)式及び(4)式のように求める方法が挙げられる。
*高吸収体領域を透過しなかった部分では
F(P1) =f(P) ・・・(3)
*高吸収体領域を透過した部分Lでは
F(P1) =α×F(G[f(P)]) ・・・(4)
即ち、高吸収体領域を透過した部分Lに関して修正用の順投影データF(G[f(P)])を所定倍(α倍)したものに置換する方法である。
【0009】
上記の方法により、修正投影データF(P1) では、高吸収体領域を透過した部分Lでは実測投影データf(P)の画素値の替わりに、修正用の順投影データF(G[f(P)])の画素値或いは高吸収体領域を透過した修正用の順投影データF(G[f(P)])に応じた画素値が利用される。従って、高吸収体領域を透過しなかった部分では実測投影データf(P)の画素値のままで、高吸収体領域を透過した部分Lのみをフィルタリング及び逆投影(FBP法による再構成)処理を行いさらに順投影したデータを置換していることになる。ここで、高吸収体領域を透過した部分Lでは実測投影データf(P)よりも、修正用の順投影データF(G[f(P)])の方が、データとして正確である。従って、修正投影データF(P1) では、アーティファクト部分が低減するような修正用の順投影データF(G[f(P)])を利用することになる。しかしながら、IRR法では修正用のデータが正確であることを要求するので、効果が限定的で、低減効果が乏しいという難点がある。一方で、アーティファクトの低減効果がより高い方法として、下記のART(Algebraic Reconstruction Technique)/EM(Expectation and Maximization)法がある。
【0010】
続いて、ART/EM法について、図10と図11とを参照して、説明をする。なお、図10は推定画像b0を1回更新して推定画像b1を求めるまでの方法を示した図であり、図11は複数回(ここではk回)更新された推定画像bk をさらに1回更新して推定画像bk+1 を求めるまでの方法を示した図である。なお、図10中のb0は任意の値で初期化される推定画像であり、b1は推定画像b0を1回更新したときの推定画像である。また、図11中のbk は、k回更新したときの推定画像であり、bk+1 は、(k+1)回更新したときの推定画像である。
【0011】
図10に示すように、上記のART/EM法の具体例では、推定画像b0を設定する。この推定画像b0はオペレータが任意の値で初期化設定することができて、例えば、全領域のX線の画素値を一様にしても構わない。そしてオペレータは断層平面G上の実測断層画像G[f(P)] を参照しながら、高吸収体に相当する部分を高吸収体領域Mとして断層平面G上に設定する(図10中の▲1▼)。そして、高吸収体領域を透過した部分Lが投影平面F上に作成される。推定画像b0を順投影して、投影平面F上に推定投影データF(b0) を得る。そして、得られた推定投影データF(b0) と実測投影データf(P)とを比較する(図10中の▲2▼)。その比較の際、高吸収体領域を透過した部分Lについては比較を行わずに無視する。つまり高吸収体領域を透過しなかった部分についてのみ、推定投影データF(b0) と実測投影データf(P)との比較を行う。これは、後述する更新についても同様に、高吸収体領域を透過しなかった部分についてのみ、推定画像の更新を行なう。
【0012】
なお、ART法は実測投影データf(P)と推定投影データF(b0) とを比較することにより、両データの差分を逆投影することに基づいて得られる比較参照画像を求める。そしてその比較参照画像によって、推定画像b0の更新を行う。一方、EM法は推定投影データF(b0) に対する実測投影データf(P)の比を逆投影することに基づいて得られる比較参照画像を求める。そして同様にその比較参照画像によって、推定画像b0の更新を行う。上記の方法によって1回更新された推定画像は推定画像b1となる(図10中の▲3▼)。また、高吸収体領域を透過した部分Lに関する画素値のデータについては、先述したように投影データ間の比較及び更新を行わずに無視するので、未定義であり区別はされない。
【0013】
推定画像b1についても、推定画像b0の更新方法と同様に、実測投影データf(P)と比較を行い、推定画像b1の更新を行う。以下、同様に推定画像を複数回更新して目的の断層画像を得る。即ち、推定画像bk の更新時は、図11に示すように、実測投影データf(P)と比較を行って推定画像bk+1 を導出する。
【0014】
上記の方法によって得られた断層画像は、IRR法と比べて、アーティファクト部分の低減効果はより高い。特に、ストリークアーティファクトについては、ART/EM法では低減効果は顕著である。しかしながら、ART/EM法では、先述したように投影データ間の比較及び更新時に高吸収体領域を透過した部分Lを無視するので、高吸収体に関するデータは未定義である。従って、高吸収体自身及びその周辺部の形状が再現できないという難点がある。また推定画像の更新の回数、即ち反復回数が少ない場合には、低コントラストな画像しか得られず、再現性の乏しい断層画像となる。高コントラストな断層画像を得るためには、推定画像の更新の回数、即ち反復回数を増やせばよい。
【0015】
この発明は、上記の事情に鑑みなされたものであって、高吸収体の吸収または散乱によって発生した断層画像上のアーティファクトを低減させるX線CT装置の画像処理方法及びX線CT装置並びにX線CT撮影用記録媒体を提供することを課題とする。
【0016】
【課題を解決するための手段】
上記のART/EM法では、先述したように、高吸収体領域を透過した部分を無視しているので、高吸収体自身及びその周辺部の形状が再現できない。従って、高吸収体に関する部分を考慮したART/EM法が望まれる。しかし高吸収体領域を透過した部分のデータをそのまま利用すると、アーティファクトの低減効果は得られない。従って、X線が高吸収体領域を透過する長さに応じて重み付けを行うART/EM法(以下、適宜「重み付きART/EM法」と略記する)が考えられる。それでも高吸収体領域を透過した部分について考慮しているので、シェーディングアーティファクトのように複数個の高吸収体に挟まれた場合には、アーティファクトの低減効果は乏しい。
【0017】
また重み付きART/EM法では、高吸収体に関する部分について考慮した影響で、反復回数を増やすことでアーティファクトの低減効果は低くなる。即ち、反復回数が少ないと、アーティファクト部分は低減されるが、低コントラストな画像しか得られない。また反復回数が多いと、高コントラストな画像が得られるが、アーティファクトの低減効果は低くなる。しかしながら、IRR法によって、高吸収体領域を透過した部分に関して実測投影データを、反復回数が少ない場合において重み付きART/EM法で得られたアーティファクトが低減された順投影データに置換して、実測投影データの修正を行い、さらにFBP法またはART/EM法等によって再構成処理を施すと、アーティファクトが低減した高コントラストな画像が得られる。
【0018】
そこで、本発明者は、IRR法と、重み付きART/EM法とを組み合わせて、さらにFBP法またはART/EM法等によって再構成処理を施すことで、アーティファクトを低減させる断層画像を求めることに想到した。
【0019】
以上のような知見に基づいて創作された本発明は、以下のような構成を取る。
即ち、請求項1に記載の発明に係る画像処理方法は、X線CT撮影時に得られた画像データを再構成して断層画像を得る画像処理方法であって、(2)被検体を透過したX線を検出することによって得られた実測投影データにフィルタリングを施した後、逆投影して実測断層画像を再構成する第1の画像再構成過程と、(3)推定画像を任意の値で初期化して設定する推定画像設定過程と、(4)第1の画像再構成過程で導出された実測断層画像に基づき高吸収体領域を設定する高吸収体領域設定過程と、(5)推定画像を前記X線の曝射方向に順投影することによって、推定投影データを導出する推定投影データ導出過程と、(6)推定投影データと実測投影データとの差分または比を逆投影することに基づいて比較参照画像を導出する比較参照画像導出過程と、(7)推定投影データと実測投影データとの差分または比に、X線が高吸収体領域を透過するパスが長くなるに従って値が小さくなる重み付き逆投影を、導入された重み関数に基づいて行って重み付き比較参照画像を導出する重み付き比較参照画像導出過程と、(8)推定画像を重み付き比較参照画像によって更新する推定画像更新過程と、(9)前記推定画像更新過程で更新された推定画像に現れる偽像の程度に応じて、前記(5)から(8)までの各過程を1回または複数回にわたり繰り返し行って偽像が低減した推定画像に更新する繰り返し動作過程と、(10)前記繰り返し動作過程で更新された推定画像を順投影することにより更新推定投影データを導出する更新推定投影データ導出過程と、(11)高吸収体領域を透過した部分に関して実測投影データを、更新推定投影データに応じたデータに置換して実測投影データの修正を行う実測投影データ修正過程と、(12)修正された実測投影データを画像再構成して断層画像を導出する第2の画像再構成過程とを備えることを特徴とする。
【0020】
請求項2に記載の発明に係るX線CT装置は、X線CT撮影時に得られた画像データを再構成して断層画像を得るX線CT装置において、(a)被検体の周りからX線を曝射するX線曝射手段と、(b)前記X線曝射手段でX線を曝射することによって、被検体を透過したX線を検出して実測投影データを得るX線検出手段と、(c)前記X線検出手段で得られた実測投影データにフィルタリングを施した後、逆投影して実測断層画像を再構成する第1の画像再構成手段と、(d)推定画像を任意の値で初期化して設定する推定画像設定手段と、(e)第1の画像再構成過程で導出された実測断層画像に基づき高吸収体領域を設定する高吸収体領域設定手段と、(f)推定画像を前記X線の曝射方向に順投影することによって、推定投影データを導出する推定投影データ導出手段と、(g)推定投影データと実測投影データとの差分または比を逆投影することに基づいて比較参照画像を導出する比較参照画像導出手段と、(h)推定投影データと実測投影データとの差分または比に、X線が高吸収体領域を透過するパスが長くなるに従って値が小さくなる重み付き逆投影を、導入された重み関数に基づいて行って重み付き比較参照画像を導出する重み付き比較参照画像導出手段と、(i)推定画像を重み付き比較参照画像によって更新する推定画像更新手段と、(j)前記推定画像更新手段で更新された推定画像に現れる偽像の程度に応じて、前記(f)から(i)までの各手段による処理を1回または複数回にわたり繰り返し行って偽像が低減した推定画像に更新する繰り返し動作手段と、(k)前記繰り返し動作手段で更新された推定画像を順投影することにより更新推定投影データを導出する更新推定投影データ導出手段と、(l)高吸収体領域を透過した部分に関して実測投影データを、更新推定投影データに応じたデータに置換して実測投影データの修正を行う実測投影データ修正手段と、(m)修正された実測投影データを画像再構成して断層画像を導出する第2の画像再構成手段とを備えることを特徴とする。
【0021】
請求項3に記載の発明に係るX線CT撮影用記録媒体は、請求項1に記載の画像処理方法をコンピュータに実行させるためのプログラムを記録した、コンピュータが読み取り可能なX線CT撮影用記録媒体である。
【0022】
【作用】
請求項1に記載の発明の作用について説明する。
本発明に係る画像処理方法によれば、高吸収体領域を設定することによって、高吸収体領域内について考慮した推定画像を得る。実測投影データと、その推定画像をX線の曝射方向に順投影することによって得られた推定投影データとを比較して、両データの差分または比を逆投影することに基づいて比較参照画像を導出する。前記両データに対して、高吸収体領域を透過するパス(長さ)に応じた重み付け逆投影を行って重み付き比較参照画像を得る。この重み付け逆投影を行う際には、導入された重み関数を用いる。これによって、高吸収体の散乱、反射による偽像を低減した推定画像に更新する。高吸収体領域に関する更新された前記推定画像、及び更新された前記推定画像をX線の曝射方向に順投影することによって得られた更新推定投影データはアーティファクトの影響が低減されている。高吸収体領域を透過した部分に関して実測投影データを、前記更新推定投影データに応じたデータに置換して、実測投影データを修正して再構成することにより、アーティファクトが低減した高コントラストな断層画像が得られる。
【0023】
請求項2に記載の発明によれば、請求項1に記載の発明の方法が、好適に実施されるので、アーティファクトが低減した高コントラストな断層画像が得られる。
【0024】
請求項3に記載の発明によれば、請求項1に記載の発明の方法をコンピュータに実行させることによって、アーティファクトが低減した高コントラストな断層画像が得られる。
【0025】
【発明の実施の形態】
以下、実施例の説明をする前に、本発明の原理を、図1を参照して説明する。
【0026】
図1は、本発明で最終的に求められる断層画像を導出するまでの流れを概略化したものである。先ず、被検体Pの周りからX線を曝射して、被検体Pを透過したX線を検出することによって、投影平面F上に実測投影データf(P)を得る。本明細書では、これを「被検体PをX線の曝射方向に投影する」(以下、適宜「被検体Pを投影する」と略記する)ともいう。なお、従来例で先述したように、ある断層画像を数学的アルゴリズムを用いて投影データを導出することを、「断層画像を順投影する」という。被検体PをX線の曝射方向に投影することによって得られた実測投影データf(P)にFBP法による再構成処理を施して断層平面G上に実測断層画像G[f(P)] を得る。本明細書では、これを「実測投影データf(P)にフィルタリングを施して逆投影する」という。なお、実測投影データf(P)に限らず、他の投影データを上記のFBP法による再構成処理を施して断層画像を得る場合でも、本明細書では「投影データにフィルタリングを施して逆投影する」という。
【0027】
次に推定画像b0を設定する。この推定画像b0は、従来のART/EM法で先述したように、オペレータが任意の値で初期化設定することができて、例えば、全領域のX線の画素値を一様にしても構わない。また、高吸収体によるX線の吸収または散乱によって、断層平面G上にアーティファクトが生じる。このアーティファクトの原因となった高吸収体に相当する部分を高吸収体領域Mとして、オペレータは断層平面G上に設定する。そして、高吸収体領域を透過した部分Lが投影平面F上に作成される。
【0028】
図1に示すように、推定画像b0を順投影することによって、投影平面F上に推定投影データF(b0) を導出する。この実測投影データf(P)と推定投影データF(bk ) (1回も更新していないときは〔F(b0)〕)とを比較して、両データの差分または比を逆投影することに基づいて比較参照画像dk (1回も更新していないときは〔d0〕 )を導出する。なお、上記の推定投影データF(bk ) は、推定画像b0をk回更新したときの推定画像bk を順投影したものである。なお、更新方法については後で詳述する。さらに高吸収体領域を透過する長さに応じた重み付き逆投影を行って重み付き比較参照画像ek (1回も更新していないときは〔e0〕 )を導出する。推定画像bkは重み付き比較参照画像ek (1回も更新していないときは〔e0〕 )によって更新されて、推定画像bk+1 (推定画像b0の更新のときはb1)を得る。即ち、これが更新である。推定投影データのF(b0) の導出から推定画像bk+1 に更新するまでが、重み付きART/EM法であり、本発明における(5)の推定投影データ導出過程から(9)の繰り返し動作過程までの各過程に相当する。図1の破線部分で符号30を付した過程が重み付きART/EM法である。
【0029】
更新された推定画像bk+1 を順投影して更新推定投影データF(bk+1)を得る。高吸収体領域を透過した部分Lに関して実測投影データf(P)の画素値を、更新推定投影データF(bk+1)に応じた画素値に置換して、実測投影データf(P)の修正を行い、修正投影データF(P2) を導出する。修正投影データF(P2) の導出までが、IRR法であり、本発明における(10)の更新推定投影データ導出過程から(11)の実測投影データ修正過程までの各過程に相当する。図1の一点鎖線部分で符号31を付した過程がIRR法である。
【0030】
修正が行われた修正投影データF(P2) を、FBP法またはART/EM法等によって再構成して、断層画像P3を求める。この断層画像P3が最終的に求められる断層画像である。また断層画像P3の導出方法は、本発明における(12)の第2の画像再構成過程に相当する。なお、第2の画像再構成過程または第2の画像再構成手段によって行なわれる画像再構成法は、FBP法だけに限定されず、ART/EM法(重み付きも含む)などによる画像再構成法を含む。
【0031】
次に、本発明に係るX線CT装置の一実施例の要部構成とその動作を、図2と図3とのブロック図を参照しながら詳しく説明する。図2は、実施例に係るX線CT装置の要部構成をあらわすブロック図である。図3は、演算部とメモリ部との内部をそれぞれ具体的にあらわしたブロック図である。
【0032】
実施例のX線CT装置における、図2中のX線管1と、X線検知器2と、検出素子CHi との機能・作用については、従来技術で述べたのと同じなのでここでは割愛することにする。それ以外の構成として、実施例のX線CT装置には、演算部10(CPU:中央演算処理部も含む)を含めた制御部3と、この制御部3を操作するためのマウス、キーボード、タッチパネル等の操作部4と、データ収集手段(DAS)7と、断層画像表示部(モニタ)8と、投影平面Fや断層平面Gの座標、及び投影データや断層画像の画素値等を記憶するメモリ部20と、記録媒体25が備わっている。このX線管1は、本発明におけるX線曝射手段に相当し、X線検知器2は、X線検出手段に相当する。
【0033】
この操作部3の入力操作についてはオペレータが行う。そしてデータの読み出し、書き込み、及び設定等の命令(コマンド)を、制御部3を介して、演算部10に行わせる。また、記録媒体25には、図3のフローチャートの手順を実行させるためのプログラムが格納されていて、記録媒体25はコンピュータに内蔵されている制御部3にそのプログラムを実行させる。プログラムが実行されると、制御部3はメモリ部20に対して書き込み及び読み出しを行い、演算部10に演算処理を行わせる。なお、操作部4には、推定画像b0を任意の値で初期化設定する推定画像設定部5と、実測断層画像G[f(P)] を参照することによって高吸収体に相当する部分を高吸収体領域Mとして設定する高吸収体領域設定部6とが備わっている。この推定画像b0(の画素値)と高吸収体領域設定部6によって設定された高吸収体領域M(の座標)とは、制御部3を介して、メモリ部20に書き込まれる。この推定画像設定部5は、本発明における推定画像設定手段に相当し、高吸収体領域設定部6は、高吸収体領域設定手段に相当し、記録媒体25はX線CT撮影用記録媒体に相当する。また、本実施例では、高吸収体領域設定部6についてはオペレータが設定を行なうように構成されているが、実測断層画像G[f(P)] の画素値等に応じて自動的に高吸収体領域Mを設定するように構成されていてもよい。
【0034】
一方、データ収集手段(DAS)7は、被検体Pの周りからX線ビームBが曝射されるのに伴って、X線検知器2で検出されるX線透過データf(P)を収集する機能を備えている。データ収集手段(DAS)7によって収集されたX線透過データf(P)、即ち実測投影データf(P)は、制御部3を介してメモリ部20に書き込まれる。次に、演算部10とメモリ部20との具体的構成について、図3を参照して詳述する。
【0035】
先ず、演算部10の具体的構成について説明する。演算部10は、図3に示すように、順投影部11と、再構成部12と、置換修正部13と、比較参照画像演算部14と、重み付き演算部15と、断層画像更新部16とを備えている。
【0036】
順投影部11は、後述する断層画像メモリ部21から読み出された断層平面G上(高吸収体領域Mをも含む)の断層画像を順投影して投影平面F上(高吸収体領域を透過した部分Lをも含む)の順投影データを導出する機能を備えている。この順投影部11は、本発明における推定投影データ導出手段と更新推定投影データ導出手段とに相当する。
【0037】
再構成部12は、後述する順投影データメモリ部22から読み出された投影平面F上(高吸収体領域を透過した部分Lをも含む)の順投影データにフィルタリングを施した後、逆投影を行い再構成することによって断層平面G上(高吸収体領域Mをも含む)の断層画像を導出する機能を備えている。つまり、再構成部12は、FBP法によって断層画像を再構成する手段を意味する。この再構成部12は、本発明における第1の画像再構成手段と、第2の画像再構成手段の一部(FBP法による再構成の場合のみ)とに相当する。
【0038】
置換修正部13は、順投影データメモリ部22から読み出された高吸収体領域を透過した部分Lに関して順投影データ(の画素値)を、別の順投影データ(の画素値)或いは後者の順投影データに応じたデータ(の画素値)に置換して、前者の順投影データの修正を行う機能を備えている。この置換修正部13は、本発明における実測投影データ修正手段に相当する。
【0039】
比較参照画像演算部14は、順投影データメモリ部22から読み出された順投影データと別の順投影データとを比較して、両データの差分または比を逆投影することに基づいて得られる比較参照画像dk を導出する機能を備えている。この比較参照画像演算部14は、本発明における比較参照画像導出手段に相当する。
【0040】
重み付き演算部15は、X線が高吸収体を透過する領域の長さに応じた後述する重み付け逆投影を行う機能を備えている。詳述すると、X線が高吸収体を透過するパスが長くなるに従って、値が小さくなる重み付き逆投影を行う。この重み付き演算部15は、本発明における重み付き比較参照画像導出手段に相当する。
【0041】
断層画像更新部16は、比較参照画像メモリ部23から読み出された重み付き比較参照画像ek によって、断層画像メモリ部21から読み出された断層画像を更新する機能を備えている。この断層画像更新部16は、本発明における推定画像更新手段に相当する。
【0042】
続いて、メモリ部20の具体的構成について説明する。メモリ部20は、図3に示すように、断層画像メモリ部21と、順投影データメモリ部22と、比較参照画像メモリ部23とを備えている。
【0043】
断層画像メモリ部21は、再構成部12と断層画像更新部16とから得られた断層平面G上の断層画像の各座標や画素値と、高吸収体領域設定部6から設定された高吸収体領域Mの座標と、推定画像設定部5から設定された断層平面G上の推定画像b0の各座標や画素値とが書き込まれ、必要に応じて、即ち操作部4の読み出し命令(READコマンド)によって制御部3を介してそれらの断層平面G上の断層画像の各座標や画素値が読み出される、読み書き可能な機能を備えている。なお、読み出された断層平面G上の断層画像は断層画像表示部(モニタ)8に表示される。そして断層画像の画素値と断層画像表示部(モニタ)8の画面の色とを対応付ける。例えば、断層画像の画素値が小さいときは画面の色を薄くして、断層画像の画素値が大きいときは画面の色を濃くして、画素値に応じて画面の色の濃淡をつけて断層画像表示部(モニタ)8に表示させる。
【0044】
順投影データメモリ部22は、データ収集手段(DAS)7によって収集された投影平面F上の被検体PのX線透過データf(P)(実測投影データf(P))の各座標や画素値と、順投影部11と置換修正部13とから得られた投影平面F上の順投影データの各座標や画素値とが書き込まれ、必要に応じて投影平面F上のそれらの順投影データの各座標や画素値が読み出される、読み書き可能な機能を備えている。なお、本実施例では投影平面F上の順投影データの各座標や画素値をモニタに表示させていないが、読み出された投影平面F上の順投影データの各座標や画素値を順投影データ用のモニタに表示させても良いし、断層画像表示部(モニタ)8に重畳表示させても構わない。
【0045】
比較参照画像メモリ部23は、比較参照画像演算部14から得られた比較参照画像dk (の画素値)と、重み付き演算部15から得られた重み付き比較参照画像ek (の画素値)とが書き込まれ、必要に応じてそれらの比較参照画像(の画素値)が読み出される、読み書き可能な機能を備えている。
【0046】
続いて、以上に説明した構成を有するX線CT装置の画像処理方法並びにX線CT撮影用記録媒体において、撮影を開始してから、断層画像を求めるまでの一連の動作を、図4と図5とのフローチャートと、図6と図7との重み付きART/EM法を参照して説明する。図4は、撮影を開始してから、断層画像を求めるまでの流れを示す(前半部分)フローチャートである。図5は、撮影を開始してから、断層画像を求めるまでの流れを示す(後半部分)フローチャートである。図6は、ART/EM法及び重み付きART/EM法について、詳しく示した図である。また、図7は、重み付きART/EM法における高吸収体を透過したパスの位置と重みの関係と、パス長に関する重み関数とを示した図である。
【0047】
〔ステップS1〕データ収集手段(DAS)7によって収集された被検体PのX線透過データf(P)は、実測投影データf(P)として、制御部3を介して順投影データメモリ部22に書き込まれる。
【0048】
〔ステップS2〕順投影データメモリ部22から読み出された実測投影データf(P)は、再構成部12によって再構成される。再構成された断層画像は実測断層画像G[f(P)] として導出されて、断層画像メモリ部21に書き込まれる。なお、書き込まれた実測断層画像G[f(P)] は、制御部3を介して断層画像メモリ部21から読み出されて、断層画像表示部(モニタ)8に表示される。この実測断層画像G[f(P)] の導出及び表示は、本発明における第1の画像再構成過程に相当する。
【0049】
〔ステップS3〕マウス、キーボード、タッチパネル等の操作部4による入力操作で、オペレータが推定画像b0を任意の値で初期化設定する。推定画像の設定の際には、操作部4内の推定画像設定部5から、制御部3を介して断層画像メモリ部21に書き込まれる。この推定画像b0の設定は、本発明における推定画像設定過程に相当する。
【0050】
〔ステップS4〕操作部4による入力操作で、実測断層画像G[f(P)] に基づいてオペレータが高吸収体領域Mを設定する。この高吸収体領域Mの設定の際には、操作部4内の高吸収体領域設定部6から、制御部3を介して断層画像メモリ部21に書き込まれる。アーティファクトの原因となった所を高吸収体領域Mとしてオペレータは設定を行う。つまり、オペレータは断層画像表示部(モニタ)8に表示された断層平面G上の実測断層画像G[f(P)] を参照しながら、高吸収体に相当する部分を高吸収体領域Mとして設定する。なお、高吸収体領域Mの設定の際には、マウス、キーボード、タッチパネル等の操作部4による入力操作で断層画像表示部(モニタ)8に直接入力設定するような構成であってもよい。この高吸収体領域Mの設定は、本発明における高吸収体領域設定過程に相当する。
【0051】
〔ステップS5〕断層画像メモリ部21から読み出された推定画像b0は、順投影部11によって順投影される。順投影されたデータは推定投影データF(b0) として投影平面F上に導出されて、順投影データメモリ部22に書き込まれる。このとき高吸収体領域を透過した部分Lが投影平面F上に作成されて、高吸収体領域を透過した部分Lの各座標も順投影データメモリ部22に書き込まれる。この推定投影データF(b0) の導出は、本発明における推定投影データ導出過程に相当する。
【0052】
〔ステップS6〕比較参照画像演算部14は、順投影データメモリ部22から読み出された推定投影データF(b0) と、実測投影データf(P)とを比較することにより、両データの差分または比を逆投影することに基づいて比較参照画像d0を導出する。そして、導出された比較参照画像d0は、比較参照画像メモリ部23に書き込まれる。この比較参照画像d0の導出は、本発明における比較参照画像導出過程に相当する。続いて、比較参照画像d0の導出について、従来のART/EM法と合わせて詳述する。
【0053】
X線管1から角度θj の方向にX線ビームBを曝射したとき、図6に示すように、投影データとして投影平面F上に実測投影データf(Pj ) が得られるとともに、再構成画像として断層平面G上に画像B[f(Pj )]が作成される。この画像B[f(Pj )]は、実測投影データf(Pj ) を逆投影したものであり、フィルタリングは行なっていない。ある検出素子CHi で検出されて収集されたX線透過データの座標は、投影平面F上の座標(CHi 、θj )に対応する。同様に断層平面G上には、角度θj の方向で被検体Pの断面と検出素子CHi とを通過する軌跡線として軌跡線Hが作成される。実測投影データf(Pj ) に関する投影平面F上の座標(CHi 、θj )の画素値をf(Pij) として、推定投影データF(b)に関する投影平面F上の座標(CHi 、θj )の画素値をF(bij) とする。このとき、推定投影データF(b)に関する投影平面F上の座標(CHi 、θj )の画素値であるF(bij) の値を、実測投影データf(Pj ) に関する投影平面F上の座標(CHi 、θj )の画素値であるf(Pij) の値に近づけるべく、ART法ではf(Pij) とF(bij) との差分であるf(Pij) −F(bij) を加算値として加算が、EM法ではF(bij) に対するf(Pij) の比であるf(Pij) /F(bij) を積算値として積算が、座標(CHi 、θj )の画素値であるF(bij) に対してそれぞれ行われる。この加算値f(Pij) −F(bij) 、または積算値f(Pij) /F(bij) を逆投影して刻み数Nで割ったものが比較参照画像になる。そして更新された推定画像は、ART法では更新する前の推定画像と比較参照画像との和になり、EM法では更新する前の推定画像と比較参照画像との積になる。なお、本明細書では、CT撮影時に回転したX線管1とX線検知器2との全角度から、1回の測定毎にX線管1とX線検知器2とが回転した角度を割った値を、刻み数Nと定義付ける。例えば、CT撮影時にX線管1とX線検知器2とが回転した全角度が360°で、1回の測定がX線管1とX線検知器2とが1°刻みで行われたとき、刻み数Nは360°/1°=360になる。単位刻み数当たりにしたのは、1回の測定の回転の刻みが違うと、更新する画像の重みが変わるからである。
【0054】
また、F(bij) をf(Pij) に変えようとすると、画像では図6に示すように、断層平面G上の軌跡線Hの部分まで一斉に変わってしまう。そのため画像が収束するまで反復する必要がある。また従来のART/EM法では、先述したように、高吸収体領域を透過しなかった部分についてのみ、実測と推定との投影データとの比較及び更新が行なわれる。
【0055】
以上より、高吸収体領域を透過しなかった部分について考えた場合、比較参照画像dk は、ART法のときは、下記の(5)式のように表され、EM法のときは、下記の(6)式のように表される。
dk =1/N×B[f(P)-F(bk )] (ART法) ・・・(5)
dk =1/N×B[f(P)/F(bk )] (EM法) ・・・(6)
【0056】
従って、従来のART/EM法の場合、更新された推定画像bk+1 は、更新する前の更新推定画像bk と比較参照画像dk とを用いると、それぞれ下記の(7)式、(8)式のように求められる。
bk+1 =bk +dk (ART法) ・・・(7)
bk+1 =bk ×dk (EM法) ・・・(8)
前述したように、上記の(7)式は、更新された推定画像bk+1 は、ART法では更新する前の推定画像bk と比較参照画像dk との和になる。上記の(8)式は、更新された推定画像bk+1 は、EM法では更新する前の推定画像bk と比較参照画像dk との積になる。ただし、上記の更新された推定画像bk+1 は、あくまでも、従来のART/EM法によって更新されたものであり、後述する重み付きART/EM法によって更新された推定画像bk+1 とは、別異であることに注意されたい。
【0057】
〔ステップS7〕推定投影データF(b0) と、実測投影データf(P)とを比較することにより、両データの差分または比に対して高吸収体を透過するパスの長さに応じた重み付け逆投影を行って、重み付き比較参照画像e0を導出する。重み付け逆投影の際には、重み付き演算部15によって、重み付けされる。重み付き比較参照画像e0は、比較参照画像メモリ部23に書き込まれる。この重み付き比較参照画像e0の導出は、本発明における重み付き比較参照画像導出過程に相当する。続いて、重み付き比較参照画像e0の導出について、重み付きART/EM法と合わせて詳述する。
【0058】
アーティファクト部分の重み付けをおこなうために、図7の(b)に示すような重み関数W(L)を導入する。この重み関数W(L)は高吸収体を透過したパスの長さLに対する関数で、そのパスの位置と重みの関係については図7の(a)に示すとおりである。即ち、X線が高吸収体領域を透過しない場合、パス長をL0としたとき、L0は0となりW(L0)は1になる。もしX線が高吸収体領域を透過して、そのパス長がL1のとき、W(L1)は1よりも小さくなる。またX線が高吸収体領域を透過して、そのパス長がL2でパス長L1よりも大きいとき、即ち高吸収体の透過領域がL1よりも長いとき、W(L2)はW(L1)よりも小さくなる。以上をまとめると、0=L0<L1<L2、1=W(L0)>W(L1)>W(L2)となる。つまり高吸収体の透過領域が長くなれば、高吸収体の吸収または散乱が大きくなって、吸収または散乱によるアーティファクトは大きくなるので、重みを軽くすることを意味している。逆に、高吸収体の透過領域が短くなれば、吸収または散乱によるアーティファクトは小さくなるので、重みを重くすることを意味している。次に、重み関数W(L)(以下、適宜「W」と略記)の式を利用して重み付き比較参照画像e0と更新後の推定画像b1とを導出する場合について説明する。
【0059】
即ち、投影データがβで表されて、重み付けを考慮した逆投影(以下、適宜「重み付き逆投影」と略記)がB1(β)で表されるとき、重み関数Wと逆投影B(β)と重み付き逆投影B1(β)とは下記の(9)式の関係で表される。
B1(β)=B(Wβ)×Const. ・・・(9)
なお、Const.は定数であり、重み付き逆投影によって通常の逆投影によって求められた元の値とは違った値になってしまう。従って、元の値と同じ値にするために、定数Const.を掛けることで調整を取っている。
【0060】
以上より、高吸収体領域を透過した部分Lも含む全領域で考えた場合、重み付き比較参照画像ek は、ART法のときは、下記の(10)式のように表され、EM法のときは、下記の(11)式のように表される。
【0061】
従って、重み付きART/EM法の場合、更新された推定画像bk+1 は、更新する前の更新推定画像bk と重み付き比較参照画像ek とを用いると、それぞれ下記の(12)式、(13)式のように求められる。
bk+1 =bk +ek (ART法) ・・・(12)
bk+1 =bk ×ek (EM法) ・・・(13)
ただし、上記の更新された推定画像bk+1 は、あくまでも、重み付きART/EM法によって更新されたものであり、先述した従来のART/EM法によって更新された更新推定画像bk+1 とは、別異であることに注意されたい。
【0062】
〔ステップS8〕断層画像メモリ部21から読み出された推定画像b0は、比較参照画像メモリ部23から読み出された重み付き比較参照画像e0によって、更新される。つまり、推定画像b0は、断層画像更新部16、即ち上述した重み付きART/EM法によって更新されるということを意味する。更新された断層画像は推定画像b1として求められて、断層画像メモリ部21に書き込まれる。なお、書き込まれた更新後の推定画像b1は、制御部3を介して断層画像メモリ部21から読み出されて、断層画像表示部(モニタ)8に表示される。この更新後の推定画像b1の導出及び表示は、本発明における推定画像更新過程に相当する。
【0063】
〔ステップS9〕オペレータは断層画像表示部(モニタ)8に表示された更新後の推定画像b1を参照しながら、その更新された推定画像b1がアーティファクトを低減させる断層画像になっているか否かを判断する。ここでは前記断層画像が低コントラストな断層画像でも構わない。もし前記断層画像になっていると判断した場合には、ステップS10に跳ぶ。もし前記断層画像になっていないと判断した場合には、ステップS5に戻って推定画像b0を更新する時と同様に推定画像b1を更新するために、ステップS5からステップS9までの工程を繰り返す。以下、同様にアーティファクトを低減させる断層画像になっていると判断するまで、推定画像bk (k回更新)から、比較参照画像dk と、重み付き比較参照画像ek と、推定画像bk をさらに1回更新した推定画像bk+1 とを求めるべく、ステップS5からステップS9までの工程を繰り返す。なお、このステップS5からステップS9までの繰り返しは、本発明における繰り返し動作過程に相当する。
【0064】
上記の重み付きART/EM法で得られた推定画像bk+1 は、反復回数が少ない場合には、従来のART/EM法と比べて高吸収体自身及びその周辺部の形状について再現性があり、低コントラストであるがアーティファクトの低減効果が得られる。
【0065】
先述したように、前記推定画像bk+1 は、ステップS5からステップS9までの工程を繰り返す回数(更新する回数)、即ち反復回数が少ないと、低コントラストだがアーティファクトの低減効果が高い断層画像になる。逆に、反復回数が多いと、推定画像bk+1 は、高コントラストな画像だがアーティファクトの低減効果が低い断層画像になる。しかしながら、次のステップS10以降のIRR法を組み合わせることによって、高吸収体領域を透過した部分に関して、実測投影データf(P)を、反復回数が少ない場合において重み付きART/EM法で得られた更新推定投影データF(bk+1)に置換して、実測投影データの修正を行い、さらにFBP法またはART/EM法等によって再構成処理を施すと、アーティファクトが低減した高コントラストな画像が得られる。つまり、重み付きART/EM法のみと比べて、同程度の高コントラストな画像を得る場合において、反復回数を減らすことができるし、アーティファクト(シェーディングアーティファクトも含む)も低減できる。
【0066】
〔ステップS10〕ステップS9で求められた推定画像bk+1 は、順投影部11によって、順投影される。順投影されたデータは更新推定投影データF(bk+1)として導出されて、順投影データメモリ部22に書き込まれる。この更新推定投影データF(bk+1)の導出は、本発明における更新推定投影データ導出過程に相当する。
【0067】
〔ステップS11〕順投影データメモリ部22から読み出された実測投影データf(P)と、同じく順投影データメモリ部22から読み出された更新推定投影データF(bk+1)とを用いて、置換修正部13は、実測投影データf(P)の修正を行い、修正投影データF(P2) を導出させる。導出された修正投影データF(P2) は、再度順投影データメモリ部22に書き込まれる。詳述すると、高吸収体領域を透過した部分Lに関して前記実測投影データf(P)の画素値を、前記順投影データメモリ部22から読み出された更新推定投影データF(bk+1)に応じた画素値に置換して、実測投影データf(P)の修正を行う。この修正投影データF(P2) の導出は、本発明における実測投影データ修正過程に相当する。次に、実測投影データf(P)を修正して修正投影データF(P2) を導出する方法について、詳述を行なう。
【0068】
高吸収体領域を透過した部分Lに関して前記実測投影データf(P)の画素値を、前記順投影データメモリ部22から読み出された更新推定投影データF(bk+1)に応じた画素値に置換する例として、高吸収体領域を透過した部分Lのみの更新推定投影データF(bk+1)をそのまま埋め込む方法が挙げられる。その場合、修正が行われた修正投影データF(P2) は下記の(14)式と(15)式のように求められる。
*高吸収体領域を透過しなかった部分では
F(P2) =f(P) ・・・(14)
*高吸収体領域を透過した部分Lでは
F(P2) =F(bk+1) ・・・(15)
【0069】
また、実施例では、高吸収体領域を透過した部分Lのみの更新推定投影データF(bk+1)をそのまま埋め込んだが、高吸収体領域を透過した部分Lのみの更新推定投影データF(bk+1)の各データの値を所定倍した値に置換してもよい。例えば、0<α<1を満たすα倍のF(bk+1)の値に置換する場合、上記の(14)式は下記の(16)式のように置き換えれる。
*高吸収体領域を透過しなかった部分では
F(P2) =α×F(bk+1) ・・・(16)
上記の(16)式のように置き換えると、アーティファクトが一層低減された投影データが求められる。
【0070】
〔ステップS12〕順投影データメモリ部22から読み出された修正投影データF(P2) は、再構成部12によって再構成される。再構成された断層画像は本発明で最終的に求められる断層画像P3として導出されて、断層画像メモリ部21に書き込まれる。或いは、再度比較参照画像演算部14と重み付き演算部15と断層画像更新部16によって、修正投影データF(P2) を重み付きART/EM法による再構成処理を行なって最終的に求められる断層画像P3を導出してもよい。この断層画像P3の導出は、本発明における第2の画像再構成過程に相当する。つまり、先述したように、最終的に求められる断層画像の導出のときは、FBP法またはART/EM法等による再構成処理で求められる。
【0071】
以上のステップにより、前述の画像処理方法及びX線CT撮影用記録媒体を有するX線CT装置では、下記の様な効果がある。即ち、反復回数を抑えた重み付きART/EM法によって、高吸収体領域を透過した部分Lのデータについて考慮して、かつ重み付けを行っている。従って、低コントラストだが高吸収体及びその周辺部が再現できて、かつアーティファクトが低減された断層画像が得られる。また反復回数を抑えた重み付きART/EM法から得られた断層画像と比検体の断層画像とを用いて、IRR法は比検体の断層画像を修正している。前述の方法から、高吸収体の部分は通常のIRR法で求められた修正用のデータと比較してより正確なデータなので、効果が限定的になることはない。また、重み付きART/EM法と比べて、同程度の高コントラストな画像を得る場合において、反復回数を減らすことができるし、アーティファクト(シェーディングアーティファクトも含む)も低減できる。
【0072】
【発明の効果】
以上に詳述したように、請求項1の発明に係るX線CT装置の画像処理方法によれば、高吸収体領域を透過した部分に関して実測投影データを、重み付けを考慮した更新推定投影データに応じたデータに置換して、実測投影データを修正して再構成することにより、アーティファクトが低減した高コントラストな断層画像を得ることができる。
【0073】
請求項2の発明に係るX線CT装置によれば、請求項1に記載の発明によって、好適に実施されるので、アーティファクトが低減した高コントラストな断層画像を得ることができる。
【0074】
請求項3の発明に係るX線CT撮影用記録媒体によれば、請求項1に記載の発明の方法をコンピュータに実行させることによって、アーティファクトが低減した高コントラストな断層画像を得ることができる。
【図面の簡単な説明】
【図1】本発明で最終的に求められる断層画像を導出するまでの流れを概略化したものである。
【図2】実施例に係るX線CT装置の要部構成をあらわすブロック図である。
【図3】演算部とメモリ部との内部をそれぞれ具体的にあらわしたブロック図である。
【図4】撮影を開始してから、断層画像を求めるまでの流れを示す(前半部分)フローチャートである。
【図5】撮影を開始してから、断層画像を求めるまでの流れを示す(後半部分)フローチャートである。
【図6】ART/EM法及び重み付きART/EM法について、詳しく示した図である。
【図7】重み付きART/EM法における高吸収体を透過したパスの位置と重みの関係と、パス長に関する重み関数とを示した図である。
【図8】従来における画像処理方法の一連の流れを示した図である。
【図9】従来のIRR法における画像処理方法の一連の流れを示した図である。
【図10】従来のART/EM法における推定画像b0を1回更新して推定画像b1を求めるまでの方法を示した図である。
【図11】従来のART/EM法における複数回更新された推定画像bk をさらに1回更新して推定画像bk+1 を求めるまでの方法を示した図である。
【符号の説明】
1 …X線管
2 …X線検知器
3 …制御部
4 …操作部
5 …推定画像設定部
6 …高吸収体領域設定部
7 …データ収集手段(DAS)
8 …断層画像表示部(モニタ)
10 …演算部
11 …順投影部
12 …再構成部
13 …置換修正部
14 …比較参照画像演算部
15 …重み付き演算部
16 …断層画像更新部
20 …メモリ部
21 …断層画像メモリ部
22 …順投影データメモリ部
23 …比較参照画像メモリ部
25 …記録媒体
30 …重み付きART/EM法
31 …IRR法
B …X線ビーム
F …投影平面
G …断層平面
L …高吸収体領域を透過した部分
M …高吸収体領域
P …被検体
Z …体軸
Claims (3)
- X線CT撮影時に得られた画像データを再構成して断層画像を得る画像処理方法であって、(2)被検体を透過したX線を検出することによって得られた実測投影データにフィルタリングを施した後、逆投影して実測断層画像を再構成する第1の画像再構成過程と、(3)推定画像を任意の値で初期化して設定する推定画像設定過程と、(4)第1の画像再構成過程で導出された実測断層画像に基づき高吸収体領域を設定する高吸収体領域設定過程と、(5)推定画像を前記X線の曝射方向に順投影することによって、推定投影データを導出する推定投影データ導出過程と、(6)推定投影データと実測投影データとの差分または比を逆投影することに基づいて比較参照画像を導出する比較参照画像導出過程と、(7)推定投影データと実測投影データとの差分または比に、X線が高吸収体領域を透過するパスが長くなるに従って値が小さくなる重み付き逆投影を、導入された重み関数に基づいて行って重み付き比較参照画像を導出する重み付き比較参照画像導出過程と、(8)推定画像を重み付き比較参照画像によって更新する推定画像更新過程と、(9)前記推定画像更新過程で更新された推定画像に現れる偽像の程度に応じて、前記(5)から(8)までの各過程を1回または複数回にわたり繰り返し行って偽像が低減した推定画像に更新する繰り返し動作過程と、(10)前記繰り返し動作過程で更新された推定画像を順投影することにより更新推定投影データを導出する更新推定投影データ導出過程と、(11)高吸収体領域を透過した部分に関して実測投影データを、更新推定投影データに応じたデータに置換して実測投影データの修正を行う実測投影データ修正過程と、(12)修正された実測投影データを画像再構成して断層画像を導出する第2の画像再構成過程とを備えることを特徴とするX線CT装置の画像処理方法。
- X線CT撮影時に得られた画像データを再構成して断層画像を得るX線CT装置において、(a)被検体の周りからX線を曝射するX線曝射手段と、(b)前記X線曝射手段でX線を曝射することによって、被検体を透過したX線を検出して実測投影データを得るX線検出手段と、(c)前記X線検出手段で得られた実測投影データにフィルタリングを施した後、逆投影して実測断層画像を再構成する第1の画像再構成手段と、(d)推定画像を任意の値で初期化して設定する推定画像設定手段と、(e)第1の画像再構成過程で導出された実測断層画像に基づき高吸収体領域を設定する高吸収体領域設定手段と、(f)推定画像を前記X線の曝射方向に順投影することによって、推定投影データを導出する推定投影データ導出手段と、(g)推定投影データと実測投影データとの差分または比を逆投影することに基づいて比較参照画像を導出する比較参照画像導出手段と、(h)推定投影データと実測投影データとの差分または比に、X線が高吸収体領域を透過するパスが長くなるに従って値が小さくなる重み付き逆投影を、導入された重み関数に基づいて行って重み付き比較参照画像を導出する重み付き比較参照画像導出手段と、(i)推定画像を重み付き比較参照画像によって更新する推定画像更新手段と、(j)前記推定画像更新手段で更新された推定画像に現れる偽像の程度に応じて、前記(f)から(i)までの各手段による処理を1回または複数回にわたり繰り返し行って偽像が低減した推定画像に更新する繰り返し動作手段と、(k)前記繰り返し動作手段で更新された推定画像を順投影することにより更新推定投影データを導出する更新推定投影データ導出手段と、(l)高吸収体領域を透過した部分に関して実測投影データを、更新推定投影データに応じたデータに置換して実測投影データの修正を行う実測投影データ修正手段と、(m)修正された実測投影データを画像再構成して断層画像を導出する第2の画像再構成手段とを備えることを特徴とするX線CT装置。
- 請求項1に記載の画像処理方法をコンピュータに実行させるためのプログラムを記録した、コンピュータが読み取り可能なX線CT撮影用記録媒体。
Priority Applications (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2000106271A JP3707347B2 (ja) | 2000-04-07 | 2000-04-07 | X線ct装置の画像処理方法及びx線ct装置並びにx線ct撮影用記録媒体 |
US09/826,317 US6426988B2 (en) | 2000-04-07 | 2001-04-05 | Image processing method of X-ray CT, X-ray CT and X-ray CT image-taking recording medium |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2000106271A JP3707347B2 (ja) | 2000-04-07 | 2000-04-07 | X線ct装置の画像処理方法及びx線ct装置並びにx線ct撮影用記録媒体 |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2001286463A JP2001286463A (ja) | 2001-10-16 |
JP3707347B2 true JP3707347B2 (ja) | 2005-10-19 |
Family
ID=18619476
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2000106271A Expired - Lifetime JP3707347B2 (ja) | 2000-04-07 | 2000-04-07 | X線ct装置の画像処理方法及びx線ct装置並びにx線ct撮影用記録媒体 |
Country Status (2)
Country | Link |
---|---|
US (1) | US6426988B2 (ja) |
JP (1) | JP3707347B2 (ja) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2010100996A1 (ja) * | 2009-03-04 | 2010-09-10 | 株式会社 日立メディコ | X線ct装置及び画像修正プログラム |
Families Citing this family (52)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20040124421A1 (en) * | 2002-09-20 | 2004-07-01 | Semiconductor Energy Laboratory Co., Ltd. | Light-emitting device and manufacturing method thereof |
JP4584550B2 (ja) * | 2003-06-05 | 2010-11-24 | 株式会社日立メディコ | X線計測装置 |
ATE385333T1 (de) * | 2003-07-18 | 2008-02-15 | Koninkl Philips Electronics Nv | Metallartefaktkorrektur bei der computertomographie |
JP3838243B2 (ja) * | 2003-09-04 | 2006-10-25 | ソニー株式会社 | 画像処理方法、および画像処理装置、並びにコンピュータ・プログラム |
ATE463808T1 (de) * | 2003-10-14 | 2010-04-15 | Koninkl Philips Electronics Nv | Iterative daten-rekonstruktion |
US7254209B2 (en) * | 2003-11-17 | 2007-08-07 | General Electric Company | Iterative CT reconstruction method using multi-modal edge information |
US7023951B2 (en) * | 2003-12-09 | 2006-04-04 | General Electric Company | Method and apparatus for reduction of artifacts in computed tomography images |
DE102004001273A1 (de) * | 2004-01-08 | 2005-08-04 | "Stiftung Caesar" (Center Of Advanced European Studies And Research) | Verfahren zur Erzeugung eines Schnittbildes |
JP2007520300A (ja) * | 2004-02-05 | 2007-07-26 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | ボクセル組織クラスを分散するctにおける高減衰オブジェクトにより生じる画像全体のアーティファクトの縮小 |
US7526060B2 (en) * | 2004-03-10 | 2009-04-28 | Koninklijke Philips Electronics N.V. | Artifact correction |
WO2006040715A1 (en) * | 2004-10-11 | 2006-04-20 | Philips Intellectual Property & Standards Gmbh | Imaging system for the generation of high-quality x-ray projections |
EP1819276A1 (en) * | 2004-10-15 | 2007-08-22 | Philips Intellectual Property & Standards GmbH | Computer tomography method |
US7991242B2 (en) | 2005-05-11 | 2011-08-02 | Optosecurity Inc. | Apparatus, method and system for screening receptacles and persons, having image distortion correction functionality |
WO2006119603A1 (en) | 2005-05-11 | 2006-11-16 | Optosecurity Inc. | Method and system for screening luggage items, cargo containers or persons |
EP3067864B1 (en) * | 2005-09-26 | 2019-05-22 | Koninklijke Philips N.V. | Iterative reconstruction with enhanced noise control filtering |
US7672421B2 (en) * | 2005-10-12 | 2010-03-02 | Siemens Medical Solutions Usa, Inc. | Reduction of streak artifacts in low dose CT imaging through multi image compounding |
DE102005050917A1 (de) * | 2005-10-24 | 2007-04-26 | Siemens Ag | Verfahren und Tomographiegerät zur Rekonstruktion einer tomographischen Darstellung eines Objektes |
EP1964055A1 (en) * | 2005-12-15 | 2008-09-03 | Philips Intellectual Property & Standards GmbH | Advanced convergence for multiple iterative algorithm |
US7899232B2 (en) | 2006-05-11 | 2011-03-01 | Optosecurity Inc. | Method and apparatus for providing threat image projection (TIP) in a luggage screening system, and luggage screening system implementing same |
US8494210B2 (en) | 2007-03-30 | 2013-07-23 | Optosecurity Inc. | User interface for use in security screening providing image enhancement capabilities and apparatus for implementing same |
US8090182B2 (en) * | 2006-11-13 | 2012-01-03 | National University Corporation Kyoto Institute Of Technology | Image reconstruction device, image reconstruction method, image reconstruction program, and CT apparatus |
WO2008084352A1 (en) * | 2007-01-04 | 2008-07-17 | Koninklijke Philips Electronics N. V. | Apparatus, method and computer program for producing a corrected image of a region of interest from acquired projection data |
US8218715B2 (en) | 2007-05-31 | 2012-07-10 | General Electric Company | Methods and systems to facilitate correcting gain fluctuations in iterative image reconstruction |
US7929742B2 (en) * | 2007-09-28 | 2011-04-19 | Siemens Medical Solutions Usa, Inc. | Method and system for computed tomographic imaging |
US8135186B2 (en) * | 2008-01-25 | 2012-03-13 | Purdue Research Foundation | Method and system for image reconstruction |
JP5601675B2 (ja) * | 2008-02-29 | 2014-10-08 | ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー | X線ct装置およびプログラム |
WO2009141779A1 (en) | 2008-05-21 | 2009-11-26 | Koninklijke Philips Electronics N.V. | Imaging apparatus for generating an image of a region of interest |
DE102008045449B4 (de) * | 2008-09-02 | 2016-03-31 | Siemens Aktiengesellschaft | Verfahren zur Erstellung von computertomographischen Aufnahmen eines Patienten mit metallischen Anteilen und Computersystem zur Durchführung des Verfahrens |
JP5355036B2 (ja) * | 2008-10-28 | 2013-11-27 | 株式会社東芝 | 核医学画像撮影装置、核医学画像再構成装置及び核医学画像再構成プログラム |
JP5342228B2 (ja) * | 2008-12-25 | 2013-11-13 | 浜松ホトニクス株式会社 | 画像処理装置および3次元pet装置 |
DE102009032059A1 (de) * | 2009-07-07 | 2011-01-13 | Siemens Aktiengesellschaft | Sinogrammbearbeitung für die Metallartefaktreduktion in der Computertomographie |
US8666137B2 (en) * | 2009-09-07 | 2014-03-04 | Koninklijke Philips N.V. | Apparatus and method for processing projection data |
JP5406063B2 (ja) * | 2010-01-27 | 2014-02-05 | 株式会社日立メディコ | 再構成演算装置、再構成演算方法、及びx線ct装置 |
US9320484B2 (en) | 2011-06-30 | 2016-04-26 | Analogic Corporation | Vibration damping for apparatus comprising rotating gantry |
EP2727083B1 (en) * | 2011-06-30 | 2017-08-09 | Analogic Corporation | Iterative image reconstruction |
CN110208295A (zh) | 2011-09-07 | 2019-09-06 | 拉皮斯坎系统股份有限公司 | 整合舱单数据和成像/检测处理的x射线检查系统 |
JP5866958B2 (ja) * | 2011-10-19 | 2016-02-24 | 株式会社島津製作所 | 放射線断層撮影装置、放射線断層画像生成処理装置および放射線断層画像生成プログラム |
DE102012204019B4 (de) * | 2012-03-14 | 2018-02-08 | Siemens Healthcare Gmbh | Verfahren zur Reduzierung von Bewegungsartefakten |
WO2014033792A1 (ja) * | 2012-08-31 | 2014-03-06 | 株式会社島津製作所 | 放射線断層画像生成装置、放射線断層撮影装置および放射線断層画像生成方法 |
CN103020928B (zh) * | 2012-11-21 | 2016-06-29 | 深圳先进技术研究院 | 锥束ct系统的金属伪影校正方法 |
WO2014181383A1 (ja) * | 2013-05-10 | 2014-11-13 | 株式会社島津製作所 | 画像処理装置 |
CN103705260A (zh) * | 2013-12-26 | 2014-04-09 | 天津大学 | 一种基于光纤通信的数字医学成像设备的数据传输系统 |
KR102379067B1 (ko) * | 2014-12-01 | 2022-03-25 | 삼성전자주식회사 | 의료 영상 장치 및 의료 영상 처리 방법 |
JP6486100B2 (ja) * | 2014-12-22 | 2019-03-20 | キヤノン株式会社 | 画像処理装置、画像処理方法、およびプログラム |
JP6371467B2 (ja) * | 2015-04-01 | 2018-08-08 | 株式会社日立製作所 | X線ct装置、再構成演算装置およびx線ct画像生成方法 |
EP3338635A4 (en) * | 2015-08-17 | 2018-08-01 | Shimadzu Corporation | Image reconstruction processing method, image reconstruction processing program, and tomography device equipped with same |
EP3391341B1 (en) | 2015-12-15 | 2025-02-19 | Koninklijke Philips N.V. | Streak artifact prediction |
GB2564038B (en) | 2016-02-22 | 2021-11-10 | Rapiscan Systems Inc | Systems and methods for detecting threats and contraband in cargo |
JP6778158B2 (ja) * | 2017-07-27 | 2020-10-28 | 株式会社日立製作所 | X線ct装置、画像生成方法および画像生成プログラム |
US10789738B2 (en) * | 2017-11-03 | 2020-09-29 | The University Of Chicago | Method and apparatus to reduce artifacts in a computed-tomography (CT) image by iterative reconstruction (IR) using a cost function with a de-emphasis operator |
CN110400286A (zh) * | 2019-06-05 | 2019-11-01 | 山东科技大学 | 一种x射线ct图像中金属针的检测定位方法 |
KR102414094B1 (ko) * | 2020-08-18 | 2022-06-27 | 연세대학교 산학협력단 | Ct영상의 금속 아티팩트 보정 방법 및 장치 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6002738A (en) * | 1995-07-07 | 1999-12-14 | Silicon Graphics, Inc. | System and method of performing tomographic reconstruction and volume rendering using texture mapping |
-
2000
- 2000-04-07 JP JP2000106271A patent/JP3707347B2/ja not_active Expired - Lifetime
-
2001
- 2001-04-05 US US09/826,317 patent/US6426988B2/en not_active Expired - Fee Related
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2010100996A1 (ja) * | 2009-03-04 | 2010-09-10 | 株式会社 日立メディコ | X線ct装置及び画像修正プログラム |
JP5404767B2 (ja) * | 2009-03-04 | 2014-02-05 | 株式会社日立メディコ | X線ct装置及び画像修正プログラム |
Also Published As
Publication number | Publication date |
---|---|
JP2001286463A (ja) | 2001-10-16 |
US20010028696A1 (en) | 2001-10-11 |
US6426988B2 (en) | 2002-07-30 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP3707347B2 (ja) | X線ct装置の画像処理方法及びx線ct装置並びにx線ct撮影用記録媒体 | |
RU2529478C2 (ru) | Способ и устройство для формирования изображений в большом поле зрения, и детектирования и компенсации артефактов движения | |
JP3894993B2 (ja) | X線源に供給されるx線管電流を変調させるシステム | |
CN102918565B (zh) | 用于采用偏心平板检测器的锥形射束计算机断层摄影成像的改进重建 | |
US9036885B2 (en) | Image reconstruction in computed tomography | |
JP4820582B2 (ja) | ヘリカルマルチスライスctのための回復ノイズを伴うヘリカルウィンドミルアーチファクトを低減する方法 | |
JP5662447B2 (ja) | 関心領域画像の再構成 | |
US9498177B2 (en) | Contrast-dependent regularization strength in the iterative reconstruction of CT images | |
JP6026214B2 (ja) | 連続マルチスケール再構成において詳細画像を補うx線コンピュータ断層撮像装置(x線ct装置)、医用画像処理装置及び医用画像処理方法 | |
JP2007117740A (ja) | 対象の断層撮影画像の再構成方法および断層撮影装置 | |
JP6062250B2 (ja) | 逐次近似法を用いたx線コンピュータ断層撮影装置(x線ct装置) | |
CN102016923A (zh) | 图像伪影的减少 | |
CN105118039B (zh) | 实现锥束ct图像重建的方法及系统 | |
JP2007527253A (ja) | 冗長データを含むヘリカルコーンビーム・コンピュータ断層撮像法に対する正確な画像再構成システム及び方法 | |
CN103229212A (zh) | 具有基于恒定方差的加权因子的迭代重建算法 | |
JP4409671B2 (ja) | 断層写真画像を作る装置と方法 | |
US8335358B2 (en) | Method and system for reconstructing a medical image of an object | |
Sunnegårdh et al. | Regularized iterative weighted filtered backprojection for helical cone‐beam CT | |
Friot et al. | Iterative tomographic reconstruction with TV prior for low-dose CBCT dental imaging | |
US20150356728A1 (en) | Method and system for substantially reducing cone beam artifacts based upon adaptive scaling factor in circular computer tomography (ct) | |
JP2002034970A (ja) | マルチ・スライスct走査の螺旋再構成の方法及び装置 | |
KR20200021156A (ko) | 사전 영상을 이용하는 고품질의 4차원 콘빔 전산화 단층 촬영 시스템 | |
JP7187131B2 (ja) | 画像生成装置、x線コンピュータ断層撮影装置及び画像生成方法 | |
US20120155729A1 (en) | Methods and systems for fast iterative reconstruction using separable system models | |
EP1374178A2 (en) | Fast computer tomography method |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20041026 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20041207 |
|
A02 | Decision of refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A02 Effective date: 20050308 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20050502 |
|
A911 | Transfer to examiner for re-examination before appeal (zenchi) |
Free format text: JAPANESE INTERMEDIATE CODE: A911 Effective date: 20050523 |
|
TRDD | Decision of grant or rejection written | ||
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20050712 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20050725 |
|
R150 | Certificate of patent or registration of utility model |
Ref document number: 3707347 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080812 Year of fee payment: 3 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090812 Year of fee payment: 4 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090812 Year of fee payment: 4 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100812 Year of fee payment: 5 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110812 Year of fee payment: 6 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110812 Year of fee payment: 6 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120812 Year of fee payment: 7 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120812 Year of fee payment: 7 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130812 Year of fee payment: 8 |
|
EXPY | Cancellation because of completion of term |