TOP > 国内特許検索 > 動き追従X線CT画像処理方法および動き追従X線CT画像処理装置 > 明細書

明細書 :動き追従X線CT画像処理方法および動き追従X線CT画像処理装置

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第6044046号 (P6044046)
登録日 平成28年11月25日(2016.11.25)
発行日 平成28年12月14日(2016.12.14)
発明の名称または考案の名称 動き追従X線CT画像処理方法および動き追従X線CT画像処理装置
国際特許分類 A61B   6/03        (2006.01)
FI A61B 6/03 350U
請求項の数または発明の数 12
全頁数 30
出願番号 特願2013-550113 (P2013-550113)
出願日 平成24年12月18日(2012.12.18)
国際出願番号 PCT/JP2012/008088
国際公開番号 WO2013/094186
国際公開日 平成25年6月27日(2013.6.27)
優先権出願番号 2011276545
優先日 平成23年12月18日(2011.12.18)
優先権主張国 日本国(JP)
審査請求日 平成27年12月15日(2015.12.15)
特許権者または実用新案権者 【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
発明者または考案者 【氏名】前田 新一
【氏名】田中 匠
【氏名】石井 信
個別代理人の代理人 【識別番号】110000822、【氏名又は名称】特許業務法人グローバル知財
審査官 【審査官】福田 裕司
参考文献・文献 特開2011-156302(JP,A)
特表2005-532137(JP,A)
特表2009-500115(JP,A)
米国特許出願公開第2007/0242796(US,A1)
福田 航 他,”混合事前分布を用いたベイズX線CT”,信学技報NC2009-133,社団法人電子情報通信学会,2010年 3月,267~272頁
Danai Laksameethanasan et al.,"A Bayesian Reconstruction Method with Marginalized Uncertainty Model for Camera Motion in Microrota,IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING VOL.57,NO.7,2010年 7月,pp.1719-1728
Bing Feng et al.,"Use of Three-Dimensional Gaussian Interpolation in the Projector/Backprojector Pair of Iterative Re,IEEE TRANSACTIONS ON MEDICAL IMAGING VOL.25,No.7,2006年 7月,pp.838-844
調査した分野 A61B 6/03
JSTPlus/JMEDPlus/JST7580(JDreamIII)
IEEE Xplore
特許請求の範囲 【請求項1】
被計測対象物X線吸収係数に関する事前知識を設定し統計推定を行うX線CT画像処理方法であって、
被計測対象物が時間と共に滑らかに変化するとし前記事前知識を第1の確率分布で設定するステップ
前記被計測対象物のX線吸収係数に対してどのようなX線投影像が観測されるのが尤もらしいかを第2の確率分布で設定するステップと、
上記の第1の確率分布と第2の確率分布の両方に依存した統計推定を行うステップと、
を備えることを特徴とする動き追従X線CT画像処理方法。
【請求項2】
被計測対象物を構成する複数組織の特性に基づいて表現されるX線吸収係数に関する事前知識を設定し統計推定を行うX線CT画像処理方法であって、
前記特性が動きの連続性を含み、
時間と共に変化する前記被計測対象物に関する前記事前知識を第1の確率分布で設定するステップと、
前記被計測対象物のX線吸収係数に対してどのようなX線投影像が観測されるのが尤もらしいかを第2の確率分布で設定するステップと、
上記の第1の確率分布と第2の確率分布の両方に依存した統計推定を行うステップと、
を備えることを特徴とするX線CT画像処理方法。
【請求項3】
1)前記第1の確率分布が、
再構成画像の画素毎に定義されるパラメータであって、
再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる確率分布で表現され、
2)前記第2の確率分布が、
再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した確率分布で表現される、
ことを特徴とする請求項1又は2に記載の動き追従X線CT画像処理方法。
【請求項4】
前記第1の確率分布は、再構成画像の隣接画素の空間的相関によって更に表現されたことを特徴とする請求項に記載の動き追従X線CT画像処理方法。
【請求項5】
前記第1の確率分布は、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することによって更に表現されたことを特徴とする請求項に記載の動き追従X線CT画像処理方法。
【請求項6】
1)投影画像と少なくともX線強度を含む計測条件が入力されるステップと、
2)前記投影画像に観測されるフォトンに関する物理過程が確率分布で定義されるステップと、
3)被計測対象物の動きに関する第1の確率分布が、被計測対象物が時間と共に滑らかに変化するとし、その時間的に連続する度合いが定義されるパラメータによって表現される確率分布で定義されるステップと、
4)観測に関する第2の確率分布再構成画像からどのような投影像が観測されると期待されるかについて記述した確率分布で定義されるステップと、
5)上記の第1の確率分布と第2の確率分布の両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織推定が行われるステップと、
を備えたことを特徴とする動き追跡X線CT画像処理方法。
【請求項7】
被計測対象物X線吸収係数に関する事前知識を設定し統計推定を行うX線CT画像処理装置であって、
被計測対象物が時間と共に滑らかに変化するとし前記事前知識を第1の確率分布で設定する手段
前記被計測対象物のX線吸収係数に対してどのようなX線投影像が観測されるのが尤もらしいかを第2の確率分布で設定する手段と、
上記の第1の確率分布と第2の確率分布の両方に依存した統計推定を行う統計推定手段と、
を備えることを特徴とする動き追従X線CT画像処理装置。
【請求項8】
1)前記第1の確率分布が、
再構成画像の画素毎に定義されるパラメータであって、
再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる確率分布で表現され、
2)前記第2の確率分布が、
再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した確率分布で表現される、
ことを特徴とする請求項に記載の動き追従X線CT画像処理装置。
【請求項9】
前記第1の確率分布は、再構成画像の隣接画素の空間的相関によって更に表現されたことを特徴とする請求項に記載の動き追従X線CT画像処理装置。
【請求項10】
前記第1の確率分布は、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することによって更に表現されたことを特徴とする請求項に記載の動き追従X線CT画像処理装置。
【請求項11】
請求項6の動き追従X線CT画像処理方法の各ステップを、
コンピュータに実行させるための動き追従X線CT画像処理プログラム。
【請求項12】
請求項11の動き追従X線CT画像処理プログラムが搭載されたX線CT画像処理装置。
発明の詳細な説明 【技術分野】
【0001】
本発明は、被計測対象物の動きに追従できるX線CTアルゴリズム、そのアルゴリズムを有するプログラムおよび該プログラムが搭載されたX線CT装置に関する。
【背景技術】
【0002】
X線CTは、非侵襲・非破壊で内部構造を知ることができる強力な手法として医療用途、欠陥検査などの産業用途など幅広く用いられている。このX線CTによる内部構造の再構成(CT画像)は、通常、内部構造を知りたい対象物が静止していることを仮定して行われてきた。しかしながら、医療用のX線CTで対象となるもののうち、心臓・血管や肺、乳幼児などに関しては動きを完全に止めることは困難である。このような動きのある物体に対して、従来通り、静止していることを想定したX線CTアルゴリズムを適用すると、モーションアーティファクトとよばれる偽像が生じることが知られている。
日本では、ここ10年ほど、心筋梗塞、狭心症などの心疾患が死因の2位を占めており、アメリカでは心疾患が死因の1位であるなど心疾患は国内外で対応が求められている疾患であり、その最も有効な対応策の一つが心疾患の異常の早期発見にあるとされるため、特に心臓のX線CTにおける正確なCT画像再構成が求められている。
【0003】
従来から、モーションアーティファクトを低減させるための技術は、その応用上の重要性から主に心臓を対象として提案されている。例えば、図14に示すように、心臓を対象としたX線CTにおいて、撮像対象を高速にスキャンすることによって撮像時間を短くし、モーションアーティファクトを低減する先行技術1(例えば非特許文献1を参照)が知られている。また、図15に示すように、心電図に基づく心位相の同期によりCT画像を再構築する先行技術2(例えば非特許文献2~4を参照)が知られている。さらに、心臓の動き推定に基づく動き修正することによりCT画像を再構築する先行技術3(例えば非特許文献5~9を参照)が知られている。
【0004】
一方、X線は強力であれば、再構成画像のSN比が高いことから、X線被曝量の低減とSN比はトレードオフの関係がある。より低いX線被曝量で従来と同程度のSN比の再構成画像を得る技術とは、従来と同じX線強度,撮像枚数でも解像度の高い画像が得られるという解像度向上を図る技術や、従来と同じX線強度や解像度でも、より撮像枚数を減らして撮像時間を短縮できる撮像時間の短縮を図る技術や、従来と同じ解像度でもX線強度を弱めて被曝量を減らせる被曝量の低減を図る技術、従来と同じX線強度や解像度でも再構成画像に含まれるノイズを低減できるノイズ除去を図る技術、従来と同じX線強度や解像度でも再構成画像のコントラスト分解能(再構成画像を何階調で表せるか、すなわち、X線吸収係数の精度)を向上するコントラスト分解能の向上を図る技術がある。
【0005】
また、上記の従来技術では、統計推定を行うものと統計推定を行わないものに大きく2つに分類できる。統計推定を行わない従来技術としては、画像再構成法の主流であるフィルター補正逆投影法(FBP:Filtered Back Projection)が知られている。FBPは、撮像対象に対して色んな方向からX線を照射して、その投影像(X線の吸収像、投影像のセットをシノグラムと呼ぶ)を得た後、得られた投影像を逆方向に投射して画像再構成する際に、シノグラムに周波数領域上でフィルターをかけボケ除去を行う画像再構成法である。
【0006】
一方、CT画像のように撮像対象の物体を異なる位置・角度から撮影した投影画像を複数枚用意できる場合、それら複数の投影画像の情報を処理することで、元の物体の断面像を再構成できることが知られている。投影画像を使用して元の物体の断面像を復元する画像再構成法のうち統計推定を行なう従来技術としては、MAP(maximum a posteriori)推定法,ベイズ推定法が知られている。
【先行技術文献】
【0007】

【非特許文献1】H. I. Goldberg etal., “Evaluation of ultrafast CT scanning of the adult abdomen,” Invest. Radiol., 24,537-543, 1989.
【非特許文献2】C. C. Morehouse et al., ‘‘Gated cardiac computed tomography with amotion phantom,’’ Radiology 134, 134-137, 1980.
【非特許文献3】P. M. Joseph and J. Whitley, ‘‘Experimental simulation evaluation ofECG-gated heart scans with a small number of views,’’ Med. Phys. 10, 444-449, 1983.
【非特許文献4】B. Ohnesorge et al., ‘‘Cardiac imaging by means ofelectrocardiographically gated multisection spiral CT: initial experience,’’ Radiology217, 564-571 , 2000.
【非特許文献5】C. J. Ritchie et al.,‘‘Correction of computed tomography motionartifacts using pixel specific back-projection,’’ IEEE Trans. Med. Imaging 15,333-342, 1996.
【非特許文献6】G. Wang et al., “A knowledge-based cone-beam x-ray CT algorithm fordynamic volumetric cardiac imaging,” Med Phys. 29(8), 1807-1822, 2002.
【非特許文献7】K. Taguchi et al., “Toward time resolved cardiac CT images withpatient dose reduction: image-based motion estimation”, Nuclear ScienceSymposium Conference Record, 4, 2029-2032, 2006.
【非特許文献8】AA. Isola et al., “Fully automatic nonrigid registration-based localmotion estimation for motion-corrected iterative cardiac CT reconstruction,” MedPhys. 37(3), 1093-1109 2010.
【非特許文献9】AA. Isola et al., “Motion compensated iterative reconstruction of aregion of interest in cardiac cone-beam CT” Comput. Med. Imag. Grap. 34(2), 149-1592010.
【発明の概要】
【発明が解決しようとする課題】
【0008】
上記の先行技術1のように、X線源を高速に回転させる場合、遠心力が大きくなり物理的な力に耐えるように設計することが困難であり、その回転速度の高速化は頭打ちとなってきているのが現状である。現在は半周で約0.25秒(1周約0.5秒)の速度が達成できているものの、心拍は約1秒に1周期であり、約0.25秒の時間に心臓の体積は大きく変わることから、先行技術1のように、心臓部分が動かないと想定するにはどうしても無理が生じるという問題がある。
【0009】
また、上記の先行技術2のように、心臓のX線CTを行う際に心電図を同時に撮り、心拍が同位相となるタイミングの投影像のみを用いて特定の位相のみの画像再構成を行う技術の場合、特定のタイミングにおける投影像のみを用いるために、結果的に被曝線量が通常よりも多くなってしまうといった問題がある。
【0010】
また、上記の先行技術3のように、心電の撮れるタイミングで同一の心位相の投影像を取得し、取得投影像から画像再構成したX線CT画像をもとに、他の心位相でのCT画像を動きがあるものとして推定する技術の場合、心臓の一周期の長さは変わることを許容するものの、一周期内では毎回、同じ体積変化を経るものと仮定して任意の時刻における心位相を特定する。しかしながら、先行技術3では、実際には動きの早い部分などでは正確に位相を合わせることが困難であり、そもそも心臓が毎回、同じようには動いていないなどの問題がある。また、血管等は心拍ごとに位置ズレが生じる場合もあり、特に心拍が短いときと長い時とでは心臓の体積変化は異なる経過を経ることが知られており、先行技術3では、これらの要因によるアーティファクトが生じるという問題がある。また、呼吸を止めることのできない患者に対しては呼吸による動きが生じてしまうためモーションアーティファクトを防げないという問題も存在する。
【0011】
そこで、本発明は、先行技術を用いた画像再構成と比べて、動きとCT画像を同時に推定して、モーションアーティファクトを低減できる、動き追従X線CT画像処理方法および動き追従X線CT画像処理装置を提供することを目的とする。
【課題を解決するための手段】
【0012】
上記状況に鑑みて、本発明の動き追従X線CT画像処理方法は、被計測対象物の動きとX線吸収係数に関する事前知識をおいて統計推定を行うX線CT画像処理方法であって、被計測対象物が時間と共に滑らかに変化することを仮定し、動きに関する第1の確率モデル(全時刻の被計測対象物の事前知識)と観測に関する第2の確率モデル(全時刻の投影像の観測モデル)を定義するステップと、上記の第1の確率モデルと第2の確率モデルの両方に依存した統計推定を行うステップと、を備える。
【0013】
かかる構成によれば、最初から統計推定を行う確率モデルとして、被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物の事前知識)と観測に関する第2の確率モデル(全時刻の投影像の観測モデル)を定義して、それらの確率モデルに基づいた推論を行なうことで、心電図などの同期を必要とせずに動きとCT画像を同時に推定することができる。
統計推定を行う確率モデルに基づくことにより、統計の分野で確立された推論手法が適用でき、またアルゴリズムの収束性が保証されるなど、優れた最適化が行いやすい。
すなわち、従来の技術においては、得られたCT再構成画像のみから動きに関するモデル(従来法は決定論的モデル)に依存して動きを推定しており、また得られた動きから観測モデルのみに依存してCT再構成画像を推定しているのに対して、本発明では、動きに関する確率的モデルおよび観測に関する確率モデルの両方を考慮して、CT再構成画像と動きの両方を同時に推定することでより精度の高い推定が行える。
また、統計推定を行う確率モデルを用いているために、他の事前知識も容易に組み込むことができ、画像再構築の品質を向上できる。
【0014】
ここで、第1の確率モデル(全時刻の被計測対象物の事前知識)は、再構成画像の画素毎に定義されるパラメータであって、再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第1の確率分布で表現され、また、第2の確率モデル(全時刻の投影像の観測モデル)は、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第2の確率分布で表現されることが好ましい。
【0015】
第1の確率モデルにおいて、再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第1の確率分布で表現されるとは、異なる時刻で異なる再構成画像が仮定され、それら異なる時刻の再構成画像が時間的に滑らかに変形して得られることが確率分布の形で表現されることである。
ここで、好ましくは、再構成画像の隣接画素の空間的相関を表現する項を加えることによって、上記の第1の確率分布をさらに特徴付ける。
また、好ましくは、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第1の確率分布をさらに特徴付ける。
【0016】
また、第2の確率モデルにおいて、各時刻の再構成画像が仮に既知であると仮定した時にその再構成画像からどのような投影像が観測されると期待されるかについての情報が確率分布の形で表現される。
【0017】
また、本発明の動き追従X線CT画像処理方法において、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第1の確率分布をさらに特徴付けることが好ましい。
さらに好ましくは、本発明の動き追従X線CT画像処理方法において、再構成画像の画素毎に与えられる隠れ状態変数によって、組織クラスを表現し、上記の第1の確率分布をさらに特徴付ける。すなわち、再構成画像をさらに特徴付けるために組織クラスを表す確率変数によって、第1の確率分布をさらに特徴付ける。組織クラスを表す確率変数は、再構成画像の画素毎に、例えば、人体などの撮像対象に関し、どのような組織(筋肉などの通常細胞、脂肪などの軟細胞、骨など)がどのようなX線吸収係数を有するかを表現する。再構成画像は時間的に滑らかに変形することが第1の確率分布に表現されるが、その表現は各時刻の再構成画像が時間的に直接相関をもつことを仮定する確率分布による表現の他、各時刻の組織クラスが時間的に滑らかに変位すること仮定することで組織クラスを通じた間接的な表現を行うことができる。また、これらの組織がどの程度の割合で分布するか、また、それぞれの組織がどの程度、空間的に連続して分布しやすいか、の分布情報が確率分布の形で表現される。
なお、組織分布に関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と組織推定の精度のさらなる向上を図ることができる。
【0018】
また、本発明の動き追従X線CT画像処理方法において、画素毎に所属する組織クラスを表現する隠れ状態変数によって、第1の確率分布をさらに特徴付けることが好ましい。
画素毎に所属する組織クラスの分布を制御する隠れ状態変数により、事前に想定した組織クラスのX線吸収係数の想定値が実際の組織クラスのX線吸収係数から大きくズレたり、想定していなかった組織クラスが存在していた場合に影響を受け難くすることができる。また、この隠れ状態変数を統計推定するにより、組織依存のCT値が存在するという事前知識がどの程度、推定される解である再構成画像へ影響するか、という影響度の度合いを制御することができる。
【0019】
ここで、本発明の動き追従X線CT画像処理方法において、上記の組織クラスの事前分布としてボルツマン分布を用いることが好ましい。
各組織の事前分布としてボルツマン分布を用いることとした理由は、被計測対象物が時間と共に滑らかに変化することを仮定したことから、隣接するフレーム間でも組織は時間と共に滑らかに変化することを仮定するべきであるが、その仮定を隣接するフレーム間の組織クラスの時間相関によっての表現しやすいことに注目したものである。
その他、組織が満たすべき性質である、組織の空間的連続性や、どの組織がどの程度の割合で存在しやすいかも同時に表現できる特性をもつ。 時間相関や空間相関の強さ、標準的な存在割合を表すパラメータは、組織クラスごとに専門家による人手による調整や実際の投影データへの当てはまりの良さを基準とした学習による調整が可能である。
なお、上記の各組織のとる典型的なX線吸収係数の値を表現するパラメータ分布としてはガウス分布を好適に用いることができ、また、各組織の空間的に連続する度合いを表現するパラメータ分布としてはボルツマン分布を好適に用いることができる。組織クラスの分布を制御する隠れ状態変数を導入した際には、各組織のX線吸収係数を表現するパラメータ分布はガウス分布の重み付き和である混合ガウス分布を好適に用いることができる。
【0020】
また、上記の本発明の動き追従X線CT画像処理方法における統計推定は、具体的には、事後確率を最大化する値による推定(MAP推定)、或いは、事後確率の期待値によるベイズ推定で行う。
ここで、MAP推定は、再構成画像の各画素の領域が属する組織とX線吸収係数の両方に関する事後分布において最も尤もらしい組合せを推定する。再構成画像の各画素の領域が属する組織との組合せを推定することは、組織クラス依存性を含めることになる。この組織クラス依存性を含めることで、例えば、空気の組織クラスはピクセル同士の連結がしやすい(空気中に他の組織は入らない)などの組織クラスに依存した知識を含めることができるのである。これにより、より現実の状況に沿った分布表現が可能となり結果的に解の推定精度を高めることが期待される。MAP推定によれば、再構成画像の事後確率を求めるのに必要な積分計算が省略することができるため、最適化を速めることができ、より高速で再構成画像を求めることができる。
【0021】
また、(近似)ベイズ推定は、X線吸収係数と再構成画像の各画素の領域が属する組織の両方に関する事後分布を推定するにあたり、事後分布をよく近似可能な試験分布を用いて近似する。ベイズ推定の場合も、MAP推定の場合と同様に、組織クラス依存性を含めることができる。これにより、より柔軟性をもった分布表現が可能となる。ベイズ推定によれば、推定の不確実性を考慮に入れた推論が可能で、またパラメータの学習が容易であり、より正確にX線CT画像処理を行うことができる。
【0022】
また、本発明の他の観点によれば、本発明の動き追従X線CT画像処理方法は、下記1)~5)のステップを備える。
1)投影画像と少なくともX線強度を含む計測条件が入力されるステップ
2)投影画像に観測されるフォトンに関する物理過程が確率分布で定義されるステップ
3)被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物 の事前知識)が、被計測対象物が時間と共に滑らかに変化することを仮定して、各組織が時間的に連続する度合いが定義されるパラメータによって特徴付けられる第1の確率分布の形で定義されるステップ
4)観測に関する第2の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについての情報が第2の確率分布の形で定義されるステップ
5)上記の第1の確率モデルと第2の確率モデルの両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織クラス推定が行われるステップ
【0023】
上記2)のフォトンに関する物理過程が確率分布は、具体的には、投影画像に観測されるフォトンに関する確率モデルがポアソン分布などで表現される。これにより、より現実の観測過程に近い物理過程を表現することができ、フォトンノイズなどによる観測の不確定性を表現できる。
上記3)~5)についての説明は、上述の説明と同様である。
【0024】
次に、本発明の動き追従X線CT画像処理装置について説明する。
本発明の動き追従X線CT画像処理装置は、被計測対象物の動きとX線吸収係数に関する事前知識をおいて統計推定を行うX線CT画像処理装置であって、被計測対象物が時間と共に滑らかに変化することを仮定し、動きに関する第1の確率モデル(全時刻の被計測対象物の事前知識)と観測に関する第2の確率モデル(全時刻の投影像の観測モデル)を定義するモデル定義手段と、上記の第1の確率モデルと第2の確率モデルの両方に依存した統計推定を行う統計推定手段を備える。
【0025】
ここで、第1の確率モデル(全時刻の被計測対象物の事前知識)は、再構成画像の画素毎に定義されるパラメータであって、再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第1の確率分布で表現され、また、第2の確率モデル(全時刻の投影像の観測モデル)は、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第2の確率分布で表現されることが好ましい。
【0026】
第1の確率モデルにおいて、再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第1の確率分布で表現されるとは、異なる時刻で異なる再構成画像が仮定され、それら異なる時刻の再構成画像が時間的に滑らかに変形して得られることが確率分布の形で表現されることである。
ここで、好ましくは、再構成画像の隣接画素の空間的相関を表現する項を加えることによって、上記の第1の確率分布をさらに特徴付ける。
また、好ましくは、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第1の確率分布をさらに特徴付ける。
【0027】
また、第2の確率モデルにおいて、各時刻の再構成画像が仮に既知であると仮定した時にその再構成画像からどのような投影像が観測されると期待されるかについての情報が確率分布の形で表現される。
【0028】
また、本発明の動き追従X線CT画像処理装置において、上述の本発明の動き追従X線CT画像処理方法と同様、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第1の確率分布をさらに特徴付けることが好ましい。
さらに好ましくは、本発明の動き追従X線CT画像処理装置において、再構成画像の画素毎に与えられる隠れ状態変数によって、組織クラスを表現し、上記の第1の確率分布をさらに特徴付ける。すなわち、再構成画像をさらに特徴付けるために組織クラスを表す確率変数によって、第1の確率分布をさらに特徴付ける。組織クラスを表す確率変数は、再構成画像の画素毎に、例えば、人体などの撮像対象に関し、どのような組織(筋肉などの通常細胞、脂肪などの軟細胞、骨など)がどのようなX線吸収係数を有するかを表現する。再構成画像は時間的に滑らかに変形することが第1の確率分布に表現されるが、その表現は各時刻の再構成画像が時間的に直接相関をもつことを仮定する確率分布による表現の他、各時刻の組織クラスが時間的に滑らかに変位すること仮定することで組織クラスを通じた間接的な表現を行うことができる。また、これらの組織がどの程度の割合で分布するか、また、それぞれの組織がどの程度、空間的に連続して分布しやすいか、の分布情報が確率分布の形で表現される。
なお、組織分布に関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と組織推定の精度のさらなる向上を図ることができる。
【0029】
また、本発明の他の観点によれば、本発明の動き追従X線CT画像処理装置は、下記1)~5)の手段を備える。
1)投影画像と少なくともX線強度を含む計測条件を入力する手段
2)投影画像に観測されるフォトンに関する物理過程を確率分布で定義する手段
3)被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物 の事前知識)が、被計測対象物が時間と共に滑らかに変化することを仮定して、その時間的に連続する度合いを定義するパラメータによって特徴付けられる第1の確率分布を定義する手段
4)観測に関する第2の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第2の確率分布を定義する手段
5)上記の第1の確率モデルと第2の確率モデルの両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織クラス推定を行う手段
【0030】
上記2)のフォトンに関する物理過程が確率分布は、具体的には、投影画像に観測されるフォトンに関する確率モデルがポアソン分布などで表現される。これにより、より現実の観測過程に近い物理過程を表現することができ、フォトンノイズなどによる観測の不確定性を表現できる。
【0031】
また、本発明の動き追従X線CT画像処理プログラムは、被計測対象物の動きとX線吸収係数に関する事前知識をおいて統計推定を行うX線CT画像処理プログラムであって、コンピュータに、下記1)~5)のステップを実行させる。
1)投影画像と少なくともX線強度を含む計測条件を入力するステップ
2)前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップ
3)被計測対象物の動きに関する第1の確率モデル(全時刻の被計測対象物の事前知識)が、各時刻の再構成画像が時間的に連続する度合いが定義されるパラメータによって特徴付けられる第1の確率分布で定義するステップ
4)観測に関する第2の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時にその再構成画像から、どのような投影像が観測されると期待されるかについて記述した第2の確率分布で定義するステップ
5)上記の第1の確率モデルと第2の確率モデルの両方に依存して、事後確率の期待値推定(ベイズ推定)もしくは事後確率の最大化推定(MAP推定)により画像再構成及び組織クラス推定を行うステップ
【0032】
また、本発明は、上記のX線CT画像処理プログラムを搭載したX線CT画像処理装置を提供することができる。
【発明の効果】
【0033】
本発明の動き追従X線CT画像処理方法および動き追従X線CT画像処理装置によれば、動きとCT画像を同時に推定して、従来技術を用いた画像再構成と比べて、モーションアーティファクトを低減できるといった効果を有する。
また、被計測対象物が時間的に滑らかに変形するという仮定に基づいた推定方法であるため、心電図などによる同期を取る必要がなく、心電図などの同期が取れない予測不可能な体動や内臓運動に対応が可能な他、全時刻の投影像を利用した推定が行えるため被曝量の低減が期待される。
【図面の簡単な説明】
【0034】
【図1】通常のCTモデルと動きを想定したCTモデルの説明図
【図2】隣接フレーム間における画素の動きの説明図
【図3】本発明の動き追従X線CT画像処理方法のフロー図(1)
【図4】本発明の動き追従X線CT画像処理方法のフロー図(2)
【図5】組織クラスごとのX線吸収係数の分布図
【図6】動きと組織クラス、組織クラスの事前分布を制御する隠れ状態変数を想定したCTの確率モデルにおける確率変数間の依存関係を表す説明図
【図7】隠れ状態変数Bの導入効果の説明図
【図8】時間と共に変形する物体の説明図
【図9】時間と共に変形する物体に対するCT画像再構成の結果を示す図
【図10】心臓の動きを伴う胸部モデルの説明図
【図11】心臓の動きを伴う胸部モデルに対するCT画像再構成の結果を示す図
【図12】組織クラスの分布を制御する隠れ状態変数bの導入効果を示す図
【図13】CT画像再構成の結果を示す図
【図14】先行技術1の説明図
【図15】先行技術2の説明図

【発明を実施するための最良の形態】
【0035】
以下、本発明の実施形態について、図面を参照しながら詳細に説明していく。ただし、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更または変形が可能である。
【0036】
本発明の動き追従X線CT画像処理方法およびX線CT画像処理プログラムでは、動きとX線吸収係数に関する事前知識をおいて統計推定を用いて画像の再構成を行う。
図1を参照し、通常のCTモデルと動きを想定したCTモデルの説明を行う。図1(1)に示すように、従来のフィルター補正逆投影法(FBP)のCTモデルの場合、動きのある心臓などの被計測対象物であっても、対象は静止していると仮定して1つの対象(CT像)から投影が得られると考えていた。本発明の場合、図1(2)に示すように、被計測対象物は初めから動いているものと仮定し、投影はこれらの異なる対象から得られるものと考える。
【0037】
そこで、静止した被計測対象物、あるいは周期的に変動する被計測対象物に対して投影を行っているという仮定はおかず、被計測対象物は時間と共に滑らかに変化することを仮定し、この仮定をもとに確率モデルを構築する。
被計測対象物は時間と共に滑らかに変化することを仮定するとは、CT画像の画素毎に動き(変形)パラメータを設け、画素毎に隣接フレーム間での動き(変形)を推定することである。具体的には、図2に示すように、j番目の画素の異なる時刻間の再構成画像(k-1番目のフレームとk番目のフレーム)の間での、動きをベクトルm(k)で表す。
以下では、時間と共に滑らかに変化する被計測対象物に対して投影を行っていることを統計モデルによって記述し、この統計モデルに基づきベイズ推定を行うことで画像再構成できることを説明する。
【0038】
図3と図4に、本発明の動き追従X線CT画像処理方法の実施形態のフローを示す。
本発明の動き追従X線CT画像処理方法のフローの一実施形態としては、図3のフローに示すように、投影画像とX線強度などの計測条件を入力し(ステップS11)、投影画像の観測フォトンに関して物理モデルで表現し(ステップS12)、X線吸収係数に関する事前知識を、被計測対象物が時間と共に滑らかに変化することを仮定し、再構成画像が時間的に連続する度合いを表現するパラメータで特徴付けた第1の確率モデル(全時刻の被計測対象物の事前知識)と第2の確率モデル(全時刻の投影像の観測モデル)とで表現し(ステップS13)、事後確率の期待値推定(ベイズ推定)で、動き,組織クラスを推定および面像再構成する(ステップS14)。
【0039】
また、本発明の動き追従X線CT画像処理方法のフローの他の実施形態としては、図4のフローに示すように、画素毎に所属する組織クラスを表現する隠れ状態変数によって、上記の第1の確率分布をさらに特徴付ける(ステップS24)。
画像再構成は不良設定性を有しているため、適当な事前知識を用いて解である再構成画像(X線吸収係数)に何らかの制約を課すことで解決する。この制約のうち、本発明の肝となるのが隣接するフレーム間の再構成画像に与えられる時間的制約である。非計測対象物が時間的に滑らかに変形することを想定しているため、図2のとおり隣接するフレーム間の再構成画像も各画素jごとに動きm(k)に応じた動きを生じる。したがって、k-1番目のフレームにおける再構成画像を動きm(k)にしたがって変形させたものは、k番目のフレームにおける再構成画像と強い時間相関をもつべきである。
この再構成画像に対する時間的制約は、再構成画像間の直接の時間相関と、後述の再構成画像の背後に潜在的に存在すると仮定される物質・組織を介した間接的な時間相関の両方、もしくは一方によって確率分布の形で表現される。
【0040】
物質・組織を介した制約は以下のとおりでる。図5に示すように、人体は脂肪,水、血、筋肉,骨,金属(例えばインプラント)といった限られた物質・組織から構成され、それぞれの物質・組織(以下、単に組織とよぶ)がある典型的なX線吸収係数をとりやすいといった制約を確率分布の形で特徴づけることができる。また、被計測対象の動きに沿って、組織が変形することを隣接するフレーム間の再構成画像の時間相関と同様な形で表現することができ、この制約が確率分布の形で特徴付けられる。その他、それぞれの組織は空間的な連続性を有しており、各組織が人体に占める割合も大まかに分かっているといった知識も組織に対する制約として同様に確率分布の形で特徴付けられる。
【0041】
図6は、図1(2)を本発明の確率モデル通りに各変数間の因果関係を具体的に書いたものである。スキャン時間を十分短い時間幅のフレームに分け、各フレーム内では対象は同一のCT画像で表されるとする(図6では、これを、x(1),x(2),・・・のように表す。)。 実際のスキャンでは、0.3秒ほどの短い時間に数百から数千の投影が得られるので、各フレームのxからは複数の投影yが得られる。不良設定性に対処するため、組織クラスZを導入し、Zが与えられたもとで、X線吸収係数xはある程度決まった値を取るとする。また、Bは事前分布の形を制御する隠れ状態変数で、事前分布をよりロバストなものにしている。
【0042】
図7は、隠れ状態変数Bを導入することによる効果を示す図である。事前分布の形を制御する隠れ状態変数Bがない場合、例えば事前分布の分散を小さくして曖昧性を低くした場合は、それに依存して事前知識の与える効果も強くなってしまう。Bを導入することにより事前分布を混合分布で表すことができ、事前知識の曖昧さと効果の強さを別々に制御することができる。
【0043】
そして、事後分布に基づいた統計推定を行う。再構成画像に関する事後分布が最大となる値(MAP推定)もしくは期待値(ベイズ推定)を求めることによって画像再構成が行える。一般に事後確率の期待値を求めるベイズ推定のほうが擾乱に強い推定が行えるが、事後確率の計算には高次元の隠れ状態変数に関する積分計算が含まれるため、解析的に実行することは困難であることから、事後確率を近似的に求める近似手法を適用することでこの計算困難さを克服する。なお、事後分布推定による統計推定においては、隠れ状態変数に関する積分計算を行わず隠れ状態変数と再構成画像の両方を含んだ同時事後確率の最大化によるMAP推定を適用することでこの計算困難さを克服しても構わない。
【0044】
(問題の定式化)
先ず、X線CT画像のように、様々な方向から投影されたT個の投影データを、D={Y(1),・・・,Y(T)}と表すことにする。各々のデータY(t)は、t番目の投影によって検出器で検出されるデータの集合であり、Y(t)={y(1),・・・,y(t)}となる。y(t)は、i番目の検出器で検出された光子の数を表す。
短時間に十分、多くの投影が得られているときには、その短時間内では、X線CT画像が同一であると想定したほうが推定すべきX線CT画像を減らすことができ計算量を抑えることができる。そこで、各投影をK通りのフレームに分け、各フレーム内では同一のX線CT画像から投影像が得られると仮定する。k番目のフレームに所属する投影データのインデックスの集合をS(k)と表し、k番目のフレームにおけるX線CT画像をX(k)と表すことにする。X線は物質を透過する際に指数的に減衰することから、観測されると期待されるX線光子の平均y(t)は、下記数式1のように表すことができる。
【0045】
【数1】
JP0006044046B2_000002t.gif

【0046】
ここで、x(k)はk番目のフレームにおける撮像対象のX線吸収係数をラスタスキャンして得られるJ次元ベクトルx={x,・・・,x}のj番目のピクセルのX線吸収係数である。vは、X線源から放出される光子数(物体が何も置かれていないときに観測され得る光子数)を表す。また、lij(t)は角度θ(t)から投影された際のi番目の検出器で検出される投影線とj番目のピクセルが交差する距離に比例するものであり、lij(t)xがt番目の投影のi番目の検出器に入射するX線に対するj番目のピクセル領域の実効的なX線吸収係数を表す。
【0047】
次に、本発明の動き追従X線CT画像処理方法および動き追従X線CT画像処理プログラムにおいて、投影画像に観測されるフォトンに関する物理モデル、事前知識としての組織に関する事前分布、すなわち、撮像対象の各組織の存在割合を表現するパラメータ分布、各組織がとりやすい典型的なX線吸収係数を表現するパラメータ分布、各組織の空間的に連続する度合いを表現するパラメータ分布、各組織の時間的に連続する度合いを表現するパラメータ分布について、以下説明を行う。
【0048】
(観測されるフォトンに関する物理モデル)
X線CTにおける観測データの主なノイズの要因は、検出される光子(フォトン)のショットノイズであると考えられる。そこで、光子(フォトン)の観測データは、投影毎、検出器毎に独立なポアソン分布に従って生成されるとしてモデル化する。
なお、ポアソン分布に特に限定されるものではなく、物理過程をより良く表現できる他の物理モデルも適用可能である。
上述の数式1は、ポアソン分布の平均を表しており、数式1の結果を用いて、下記数式2で表わす。ノイズは投影、画素(ピクセル)毎に独立であることから、数式2は全体として独立なポアソン分布の積で表されている。数式2において、S(k)は、k番目のフレームに含まれる投影のインデックスの集合を表している。
【0049】
【数2】
JP0006044046B2_000003t.gif

【0050】
(動きに関する事前分布)
動きθが与えられたもとでの組織クラスの分布はボルツマン分布で表されると仮定して、組織クラスに関する事前分布を下記数式3に示すようにモデル化する。ここで、Aは正規化定数であり、H(Z)はどのような組織クラスの組み合わせが尤もらしいか表すエネルギー関数である。エネルギー関数H(Z)は、下記数式4で表される。
ここで、z(k) = [zj1(k), ・・・, zjC(k) ] ∈ {0,1}は、k番目のフレームにおいてj番目のピクセルが属する組織クラスを表すベクトルであり、組織クラスcに属する場合はzjc(k) =1、その他の成分は0となる。大文字のZは小文字の変数z(k)の集合を表す。
【0051】
【数3】
JP0006044046B2_000004t.gif

【0052】
【数4】
JP0006044046B2_000005t.gif

【0053】
上記数式4において、η(j)は、j番目のピクセルの4近傍、N(j)はj番目のピクセルを含むその周囲9×9(ピクセル)。1項目は組織の割合、2項目は同一フレームにおける組織の空間的なめらかさ(つながりやすさ)、3項目は隣接するフレーム間での組織の時間的になめらかな変位・変形を表現するものである。また、数式4中、G(i,j,m(k)(θ))はH(Z)の3項目の重みであり、下記数式5で表されるもので、隣り合うフレーム(k-1番目とk番目) において、k-1番目のフレームの各ピクセルが現在の動きm(k)(θ)にしたがって動いたとして、k番目のフレームの各ピクセルとの距離を計算し、その距離が近いほど重みが大きくなる。
【0054】
【数5】
JP0006044046B2_000006t.gif

【0055】
(再構成画像Xに関する事前分布)
再構成画像Xの事前分布を下記数6に示す。組織クラスZと隠れ状態変数Bが与えられたもとで、X線吸収係数はガウス分布として表現している。以下に説明するように、ガウス分布の平均は組織クラスごとに設定する。また、その組織クラスで典型的なX線吸収係数の値への接近の度合いは隠れ状態変数bによって制御される。各フレームk、各画素jにおいてbj1(k), bj2(k)はいずれか一方が1をとりいずれか一方が0をとる変数である。組織クラスcごとに与えられるその組織クラスcに典型的なX線吸収係数の値vc1, vc2は通常、同じ値vc1=vc2=vが設定される。たとえば、平均vの値は、各々の組織(空気,脂肪,筋肉,骨,金属)ごとにそれぞれ、v=0.000,v=0.018,v=0.022,v=0.040,v=0.120と設定する。また、分散パラメータεc1, εc2は組織クラスcごとに固定した値であるが、これは一方が大きく、一方が小さい正の値をとるパラメータである。vc1, vc2c1, εc2はいずれも専門家による人手あるいは実際に得られた多くの投影像への当てはまりの良さを基準としたパラメータ学習によって与えられる。
数式6の1項目は典型的なX線吸収係数の値への再構成画像が表すX線吸収係数の近接度の度合いを決定する。隠れ状態変数bj1(k), bj2(k)はいずれか一方が1、いずれか一方が0をとる変数であるため、再構成画像の組織クラスcの典型的なX線吸収係数への値の近接度は組織クラスcのみによって決定されるのではなく各フレームk、各画素jごとにbj1(k), bj2(k)のとる値によっても制御されることになる。このような組織クラスcごとに固定された分散パラメータεc1, εc2の他に、与えられた投影像から推定される隠れ状態変数bj1(k), bj2(k)を用意し推定することで各フレーム、各画素ごとに適応的に近接度を制御可能とし、
これによって組織クラスごとにどのX線吸収係数をとりやすいかに関する不確かさとその影響の強さを別々に制御することができる。数式6の2項目と3項目はそれぞれ再構成画像が空間的に滑らかに変化すること、動きベクトルmに基づいて時間的に滑らかに変形することを表す項であり、GspatialとGteemporalは空間相関の強さと時間相関の強さを制御するパラメータである。
【0056】
【数6】
JP0006044046B2_000007t.gif

【0057】
ここで、隠れ状態変数Bはボルツマン分布に従うと仮定した場合(下記数式7)、隠れ状態変数Bのどのような組み合わせが尤もらしいかを表すエネルギー関数H(B)は下記数式8で表される。ここで、数式8の1項目は2値の隠れ変数の集合であるBがそれぞれ2値のいずれの値をとりやすいかを表し、2項目は分布の空間的なめらかさ(隣接するピクセルでの同じ分布の取りやすさ)を制御するものである。
【0058】
【数7】
JP0006044046B2_000008t.gif

【0059】
【数8】
JP0006044046B2_000009t.gif

【0060】
また、k-1番目のフレームからk番目のフレームへのj番目のピクセルの動きを表すmj(k)(θ)は、下記数式9で表される。数式9は、少数の代表動きベクトルθi(l)(i=1,・・・,J, l=1,・・・,K)からk番目のフレームにおけるj番目のピクセルの動きベクトルmj(k)(θ)を補間により求める式である。空間的、時間的に近い代表動きベクトルθi(l)(i=1,・・・,J, l=1,・・・,K)の重みが大きくなるようになっている。なお、数式9において、rはi番目の画素の位置ベクトルを表す。
【0061】
【数9】
JP0006044046B2_000010t.gif

【0062】
(事後分布のベイズ推定)
CT画像の再構成では、X線吸収係数xと組織クラスzと隠れ状態変数bは、観測データDと被計測対象物の動きθの事後分布p(X,Z,B|D,θ)として推定される。事後分布p(X,Z,B|D,θ)は、ベイズの定理により分解でき、それぞれの事前分布と尤度の積に比例する形で、下記数式10のように表される。MAP推定においては事後分布p(X,Z,B|D,θ)の最大値を求めることにより、高次元の隠れ状態変数に関する和計算の必要がなくなる。
【0063】
【数10】
JP0006044046B2_000011t.gif

【0064】
(近似)ベイズ推定では、計算困難な事後分布p(X,Z,B|D,θ)を近似する試験分布q(X,Z,B)によって評価される。なお、ベイズ推定自体は、事後分布を計算しその期待値を推定値とする方法であって、事後分布やその期待値が解析的に計算可能なら試験分布は不要である。
試験分布q(X,Z,B)は、KL距離と呼ばれる分布の近さを表す指標を最小化するように決定される。ここでのKL距離は下記数式11で定義される。
【0065】
【数11】
JP0006044046B2_000012t.gif

【0066】
ここで、<・>q(X,Z,B)は、試験分布q(X,Z,B)に関して積分計算を行うことを表す演算子である。すなわち、試験分布q(X,Z,B)に関して期待値を取ることを意味する。KL距離は常に非負であり、q=pのときにのみ0(ゼロ)になる。事後分布の計算困難性に対処するために、計算可能な分布の集合の中において事後分布に最も近いものを探索する。
試験分布q(X,Z,B)は、KL距離(上記数式11)の最小化が可能であるようなもの、あるいは近似的な最小化が可能であるようなものであれば、任意に選ぶことができる。そのため、計算を簡単にするために、試験分布q(X,Z,B)には、下記数式12に示すように、試験分布q(X,Z,B)に対して因子化仮定と呼ばれる変数xとzとbの間の独立性の仮定を行う。すなわち、X,Z,Bは独立で、さらにフレーム、ピクセル毎にも独立であるとする。 さらに、xに関する分布q(x)がガウス分布になるものと仮定して、下記数式13に示されるようにモデル化する。これらの仮定の下で、q(X),q(Z),q(B)を交互に最適化する交互最適化を用い最適試験分布を求めている。
【0067】
【数12】
JP0006044046B2_000013t.gif

【0068】
【数13】
JP0006044046B2_000014t.gif

【0069】
q(X)とq(B)を現在、得られている推定値に固定し、q(Z)のみを上記数式11のDKLが最小化できるよう最適化する際には、解析解が得られる。同様に、q(X)とq(Z)を現在、得られている推定値に固定し、q(B)のみを上記数式11のDKLが最小化できるよう最適化する際にも、解析解が得られる。
q(Z)とq(B)を現在、得られている推定値に固定しても、q(X)は解析解が得られないので、上記数式13におけるq(X)のパラメータμ,σは数値最適化によって求める。これらの分布q(X), q(Z), q(B)は、事前に決められた所定の収束条件を満たすまで最適化が繰り返され、解が求められる。
【0070】
以上、本発明の動き追従X線CT画像処理方法および動き追従X線CT画像処理プログラムにおける処理内容について説明した。以下の実施例では、変形する物体を人工的に生成し、その投影像から画面再構成を行い、従来技術による画面再構成と比較して評価したものである。下記の評価実験から、本発明の動き追従X線CT画像処理方法および動き追従X線CT画像処理プログラムの有用性の理解が深まるであろう。
【実施例1】
【0071】
実施例1では、変形する物体を人工的に生成し、その投影像から本発明の動き追従X線CT画像処理方法により画面再構成を行う計算機シミュレーションを行った結果を説明する。
先ず、図8に示す画像は、人工的に生成した変形する物体(64×64ピクセル)であり、動きのある物体を模擬するものである。図8に示すように、それぞれの円形の画面の中央部分が物体を表しており、これが時刻と共に(a)→(b)→(c)→(d)→(e)と膨張していくものとした。物体の変形は、縦と横にそれぞれ異なる比率で膨張していくものとした。
X線CT画像の再構成の結果を図9に示す。図9は、真の画像と、動きを想定しない従来法として最もよく用いられているフィルター補正逆投影法(FBP)による画像再構成の結果と、本発明による画像再構成(提案法)の結果とを比較したものである。
それぞれの手法による再構成画像と真の画像との差をPSNR(Peak signal to noise ratio)によって評価した。PSNRは下記数式14のように定義されるものである。
【0072】
【数14】
JP0006044046B2_000015t.gif

【0073】
上記数式14において、MSEは真の画像と再構成画像の画素値間の平均二乗誤差を表し、maxは、真の画像の画素値の最大値を表す。図9の各図の下部にPSNRの値を付記している。
図9の結果から、本発明はPSNRで平均して約3.8(dB)の改善をもたらしていることがわかる。約3.8(dB)の改善は平均二乗誤差でいうと約42%の誤差を低減していることを意味する。得られる画像からも中心の変形する物体の動きをよく表現できていることがわかる。
【0074】
次に、図10に示す画像は、心臓の動きを伴う胸部モデルの画像であり、心臓の動きを伴う胸部を模擬するものである。図10に示すように、心臓に対応する中心の上の楕円がだんだん大きくなるように、すなわち時刻と共に(a)→(b)→(c)→(d)→(e)と心臓に対応する楕円が膨張していくものとした。楕円の変形は、縦と横にそれぞれ異なる比率で膨張していくものとした。胸部モデルの画像サイズは、256×256(ピクセル)、投影数は360°の間に360枚とした。
X線CT画像の再構成の結果を図11に示す。図11は、真の画像と、動きを想定しない従来法として最もよく用いられているフィルター補正逆投影法(FBP)による画像再構成の結果と、本発明による画像再構成(提案法)の結果とを比較したものである。
図11の結果から、本発明はPSNRで平均して約1.1(dB)の改善をもたらしていることがわかる。約1.1(dB)の改善は平均二乗誤差でいうと約22%の誤差を低減していることを意味する。得られる画像からも心臓の動きをよく表現できていることがわかる。
【実施例2】
【0075】
本発明の動き追従X線CT画像処理方法では、上述したように、統計推定を行う確率モデルを用いているために、他の事前知識も容易に組み込むことができ、画像再構築の品質を向上できる。例えば、各組織の代表的なX線吸収係数を考慮に入れて推定を行うことで推定精度を向上させることができる。この際、各組織クラスの分布を制御する隠れ状態変数を含めることで推定精度の向上、ロバストな推定ができる。その効果を調べた実験について、以下説明する。
【0076】
図12は、想定した組織クラスのX線吸収係数と、実際に実験で用いた物体のもつX線吸収係数の値を示している。ここでは、空気、脂肪、筋肉、骨の4つの組織クラスが存在すると仮定し、それぞれの組織クラスで代表的なX線吸収係数があるものとしている。
実際の投影では、想定した代表値からX線吸収係数の値のズレが生じることを考えて、図中に符号aで示した線のように、脂肪と骨のクラスのX線吸収係数は10%だけX線吸収係数が高いものとした。また図中に符号bで示した線は、その他の外れ値のとるX線吸収係数の値である。各組織クラスのうち、脂肪と筋肉に関して想定した代表的なX線吸収係数から5%大きな値をとる物体に対してX線CTを行った。
【0077】
図13は、上記のような条件設定で、X線CT再構成の結果を比較したものである。図13において、(a)は真の画像であり、(b)はFBPによる再構成画像であり、(c)は本発明による再構成画像で隠れ状態変数Bの導入が無しのものであり、(d)は本発明による再構成画像で隠れ状態変数Bの導入が有りのものである。また(e)は、推定された隠れ状態変数Bを表している。
隠れ状態変数Bを導入した場合と導入しない場合とで、PSNRで0.8(dB)の改善(平均二乗誤差で17%程度の改善)となった。従来法であるフィルター補正逆投影法(FBP)と比較すると17.6(dB)の改善となっている。また、各組織におけるX線吸収係数が想定からズレた場合を考えているが、5%と小さいシフトしか想定していない。このため、より大きなズレが生じる場合には、より大きな効果が得られると期待される。
【0078】
以上の結果は、比較的、確率モデルを簡単化したものでの結果であり、推定精度には向上の余地が存在するものの、推定精度自体は実用化のための十分なレベルを達成していると考える。なお、実用化のための最適なパラメータ設定は、実際の計測機器の仕様や測定対象の性質に依存して決定されるべきものであるため、これらの決定のためには実際の計測機器の仕様や測定対象の性質を知る必要がある。
【実施例3】
【0079】
次に、実施例3では、上記実施例1,2と異なる動き追従X線CTアルゴリズムについて説明する。
(1)確率モデルについて
(1-a)尤度
再構成の対象となる物体のCT画像を1次元ベクトルX=[x,・・・,x](ここで,Jは再構成画像の画素数)で表現する。 ここでは、時間に応じて変形もしくは移動する物体に対する投影を考えるため、時間をK個のフレームS(k), k∈{1,・・・,K}に分割し、k番目のフレームでのCT画像をX(k)=[x(k),・・・,x(k)]で表現する。全フレームのX(k)(k∈{1,・・・,K})を、纏めてX={X(1),・・・,X(k)}で表現する。
時刻tにおける物体に対する投影による像をY(t)=[y(t),・・・,y(t)] (ここで、Iは検出器数)とし、全時刻の投影を、纏めてD={Y(1),・・・,Y(T)}で表現する。時刻tがt∈S(k)を満たす場合、CT画像X(k)の投影を得たと想定するので、光電変換で生じるノイズが独立なポアソン分布に従うことを考えたとき、下記数式15,数式16で表現できる。
【0080】
【数15】
JP0006044046B2_000016t.gif

【0081】
【数16】
JP0006044046B2_000017t.gif

【0082】
上記数式(16)におけるL(t)は、mxnの行列であり、その(i,j)要素は照射されたX線のうち検出器iに入射する光線が、再構成画像の画素jとどれだけ交差するかを表している。また、上記数式(16)におけるbは、ブランクスキャン時のフォトン数(照射されたフォトン数)である。ここで、ノイズの分布にはガウス分布を採用することも可能である。
【0083】
(1-b)事前分布
事前分布を下記数式17で表現する。
【0084】
【数17】
JP0006044046B2_000018t.gif

【0085】
ここで、Z,Wは各フレームのCT画像の画素それぞれに対応する変数からなり、それぞれ組織クラスを表す変数Z={Z(1),・・・,Z(K)},
分布形状を制御する変数B={B(1),・・・,B(K)}である。さらに、k番目のフレームにおける組織クラス, 分布形状を制御する変数は、それぞれZ(k)=[z(k),・・・,z(k)],B(k)=[b(k),・・・,b(k)]と表現される。想定する組織クラスの数をCとしたとき、各画素の組織クラスは、C次元ベクトルz(k)=[zj1(k),・・・,zjC(k)]で表現される。zjC(k)は、そのk番目のフレームの画素jがとる組織クラスcに対してのみzjC(k)=1をとり、その他の要素はゼロをとる。
したがって、下記数式18の制約を満足する。また、k番目のフレームの画素jの分布形状を制御する変数b(k)は、は、下記数式19の2次元ベクトルとしている。なお、より複雑な分布を表現するためには、dの取りうる値の数を2より大きな任意の数に増やすことで対応可能である。
【0086】
【数18】
JP0006044046B2_000019t.gif

【0087】
【数19】
JP0006044046B2_000020t.gif

【0088】
パラメータθは、動きを表すパラメータとして定義され、事前分布は、下記数式20で表現される。
数式20におけるE(X,Z,B,θ)は、下記数式21で表現される。
【0089】
【数20】
JP0006044046B2_000021t.gif

【0090】
【数21】
JP0006044046B2_000022t.gif

【0091】
以下に、パラメータおよび変数の説明をする。
まず、下記数式22で表現されるパラメータは、zjc(k)=zsc´(k)=bjd(k)=bsd´(k)=1のときのみ、値Gx,cdspatialをとる。また、下記数式23で表現されるパラメータは、zjc(k)=zsc´(k)=1のときのみ、値Gx,cspatialをとる。
【0092】
【数22】
JP0006044046B2_000023t.gif

【0093】
【数23】
JP0006044046B2_000024t.gif

【0094】
したがって、Gx,cdspatialを変更することによって、隣接画素の組織クラスと分布形状制御変数が同一だった時に、画素値xがどれだけ滑らかとなるべきかを制御できる。また、Gx,cspatialを変更することによって、分布形状制御変数によらず隣接画素の組織クラスが同一だった時に画素値xがどれだけ滑らかとなるべきかを制御できる。同様にGz,cspatialは隣接画素がどれだけ同じ組織クラスとなりやすいか、Gb,cdspatialは隣接画素が同じ組織クラスだった時にどれだけ同じ分布形状制御変数をとるかを制御する。Gcdselfは組織クラス変数、分布形状制御変数のどのような組み合わせが生じやすいかを表す。Gtemporal, Gtemporalはそれぞれ変数x、zがどれだけ動きパラメータθに沿って時間的に遷移しやすいかを制御する変数である。η(j),N(j)は、画素jの近傍の画素インデックスの集合を表す。それぞれ滑らかさをどれだけの範囲に及ぼすかを決定するものであり、計算量の節約のためには小さくとることが好ましい。
【0095】
また、関数G(i,j,m(k),(θ))は、下記数式24で定義される。ここで、rは画素iの位置ベクトル、m(k)(θ)は画素iのフレームkからk+1への移動量を表す動きベクトルである。つまり、k番目のフレームの画素iは(k+1)番目のフレームにおいて、位置rから位置r+m(k)(θ)に移動すると期待されることを表現する。なお、γはある正の定数である。したがって、関数G(i,j,m(k),(θ))は、期待される移動量に近い画素ほど関数G(i,j,m(k),(θ))が大きな値をとるような重み付け関数として働くことになる。
【0096】
【数24】
JP0006044046B2_000025t.gif

【0097】
ここで、移動量を表す動きベクトルm(k)は、下記数式25で定義される。これは、exp(-α||r-r||-β|k-l|)が、時空間的な重みと考えれば、動きベクトルm(k)がパラメータθ(l)の補間で得られることを表している。α,βはそれぞれ空間的補間と時間的補間の強さを制御するパラメータである。
【0098】
【数25】
JP0006044046B2_000026t.gif

【0099】
下記数式26は、Gx,cdspatial=Gx,cspatial=Gtemporal=0が成り立つ時、下記数式27のように、P(X|Z,B,θ)が、Z,Bが与えられた下で独立という条件付き独立な分布として表現される。このとき、P(X|Z,B,θ)は数式27を改めて下記数式28のように表現できる。
【0100】
【数26】
JP0006044046B2_000027t.gif

【0101】
【数27】
JP0006044046B2_000028t.gif

【0102】
【数28】
JP0006044046B2_000029t.gif

【0103】
また、上記数式27におけるp(Z,B|θ)は下記数式29で、上記数式26におけるE(X,Z,B,θ)は下記数式30で定義される。このように定義することで設定すべきパラメータを減らすことができる。
【0104】
【数29】
JP0006044046B2_000030t.gif

【0105】
【数30】
JP0006044046B2_000031t.gif

【0106】
P(x(k)|zjc(k)=1,bjd(k)=1,θ)に関しては、CT値xが非負の値をとることを考慮して、ガンマ分布を利用することも可能である。このとき、P(x(k)|zjc(k)=1,bjd(k)=1,θ)は、下記数式31で表現される。数式31において、ucd,vcdはそれぞれガンマ分布のパラメータである。
【0107】
【数31】
JP0006044046B2_000032t.gif

【0108】
次に、上記の確率モデルを用いて、動きのパラメータθを推定しつつ、CT画像X, 組織クラス変数Z, 分布形状制御変数Bを推定する方法について説明する。
(2)自由エネルギー最小化に基づく推定方法について
表記を簡単にするために観測変数yから隠れ変数xを推定する一般的な問題を考える。上記の確率モデルでは、観測変数をD,隠れ変数をX,Z,Bと考えれば対応がつく。また、自由エネルギーFは、下記数式32で定義される。数式32において、q(x)は試験分布と呼ばれる分布である。
【0109】
【数32】
JP0006044046B2_000033t.gif

【0110】
対数尤度logp(y|θ)は、自由エネルギーFを用いると下記数式33で表現できる。ここで、KL[q(x)|p(x)]は、xが実数の場合は数式34で定義される関数であり、xが離散変数の場合は数式35で定義される関数であり、確率分布間の近さを測る尺度として用いられる。これは常に非負の値をとり、q(x)=p(x)が成り立つときのみゼロをとる。
【0111】
【数33】
JP0006044046B2_000034t.gif

【0112】
【数34】
JP0006044046B2_000035t.gif

【0113】
【数35】
JP0006044046B2_000036t.gif

【0114】
また、パラメータ推定では、しばしば最尤推定が用いられる。最尤推定の目的は下記数式36の最適化である。上記の式変形からわかるように対数尤度は自由エネルギーを用いて、下記数式37と表現できる。ここで、数式37における左辺は、試験分布q(x)に依存しない量であるので、右辺も試験分布のとり方に依存しない。
【0115】
【数36】
JP0006044046B2_000037t.gif

【0116】
【数37】
JP0006044046B2_000038t.gif

【0117】
また、数式38であるから、数式39が成り立つ。この式から右辺は、KL[q(x)|p(x|y,θ)]=0のとき、すなわち、q(x)=p(x|y,θ)のとき、最小値-logp(y|θ)が達成されることがわかる。つまり、対数尤度を最大化するパラメータは下記数式40と表現できる。
【0118】
【数38】
JP0006044046B2_000039t.gif

【0119】
【数39】
JP0006044046B2_000040t.gif

【0120】
【数40】
JP0006044046B2_000041t.gif

【0121】
したがって、左辺の自由エネルギーを試験分布q(x)に関して最小化したうえで、パラメータθに関して最小化することで最尤推定が行えることがわかる。この事実より、対数尤度の最大化のかわりに自由エネルギーを試験分布q(x)とパラメータθについて最大化するという方法が考えられる。
ここでは、数式41を求めるために、数式42に示す2つのステップを繰り返す。
試験分布q(x)の最適値は、事後分布p(x|y,θ)で与えられるが、一般に事後分布の解析的な計算は困難である。そこで、q(x)の候補を制約することで最適化を可能にする。
【0122】
【数41】
JP0006044046B2_000042t.gif

【0123】
【数42】
JP0006044046B2_000043t.gif

【0124】
ここでは、下記数式43なる独立性を仮定した。この仮定のもとで、交互最適化によって最適化することを考えると、下記数式44になる。
そして、変分法によって、数式45の解を解析的に求めることができる。一方、数式46に関しては解析的に求まらないので、q(x(k))として、下記数式47のガウス分布を仮定するか、或いは、下記数式48のガンマ分布を仮定する。
【0125】
【数43】
JP0006044046B2_000044t.gif

【0126】
【数44】
JP0006044046B2_000045t.gif

【0127】
【数45】
JP0006044046B2_000046t.gif

【0128】
【数46】
JP0006044046B2_000047t.gif

【0129】
【数47】
JP0006044046B2_000048t.gif

【0130】
【数48】
JP0006044046B2_000049t.gif

【0131】
上記のようにすることで、数式49を解析的に評価することができる。また、数式49において、q(X)のパラメータ(q(x(k))が、ガウス分布に従うとすれば平均μと分散σ、ガンマ分布に従うとすればパラメータ(ρcd,ωcd)に関して共役勾配法などの非線形最適化法で最適化することができる。なお、ガンマ分布はxの非負性を表現できる好ましい性質を有する。
【0132】
【数49】
JP0006044046B2_000050t.gif

【産業上の利用可能性】
【0133】
本発明は、医用ならびに産業用のX線CT装置に有用である。医用画像とくに心臓のX線CTが重要な応用であるが、産業用の非破壊検査においても微小な物体の測定時には、ステージの微小な動きが問題になったり、流体が存在する場合にはその流体の変形が問題になったりすることから、様々な用途において、本発明のモーションアーティファクトの低減化が役立つことが期待される。

図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14