(57)【特許請求の範囲】
【請求項１】
被計測対象物のＸ線吸収係数に関する事前知識を設定し統計推定を行うＸ線ＣＴ画像処理方法であって、
被計測対象物が時間と共に滑らかに変化するとし、前記事前知識を第１の確率分布で設定するステップと、
前記被計測対象物のＸ線吸収係数に対してどのようなＸ線投影像が観測されるのが尤もらしいかを第２の確率分布で設定するステップと、
上記の第１の確率分布と第２の確率分布の両方に依存した統計推定を行うステップと、
を備えることを特徴とする動き追従Ｘ線ＣＴ画像処理方法。
【請求項２】
被計測対象物を構成する複数組織の特性に基づいて表現されるＸ線吸収係数に関する事前知識を設定し統計推定を行うＸ線ＣＴ画像処理方法であって、
前記特性が動きの連続性を含み、
時間と共に変化する前記被計測対象物に関する前記事前知識を第１の確率分布で設定するステップと、
前記被計測対象物のＸ線吸収係数に対してどのようなＸ線投影像が観測されるのが尤もらしいかを第２の確率分布で設定するステップと、
上記の第１の確率分布と第２の確率分布の両方に依存した統計推定を行うステップと、
を備えることを特徴とするＸ線ＣＴ画像処理方法。
【請求項３】
１）前記第１の確率分布が、
再構成画像の画素毎に定義されるパラメータであって、
再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる確率分布で表現され、
２）前記第２の確率分布が、
再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した確率分布で表現される、
ことを特徴とする請求項１又は２に記載の動き追従Ｘ線ＣＴ画像処理方法。
【請求項４】
前記第１の確率分布は、再構成画像の隣接画素の空間的相関によって更に表現されたことを特徴とする請求項３に記載の動き追従Ｘ線ＣＴ画像処理方法。
【請求項５】
前記第１の確率分布は、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することによって更に表現されたことを特徴とする請求項３に記載の動き追従Ｘ線ＣＴ画像処理方法。
【請求項６】
１）投影画像と少なくともＸ線強度を含む計測条件が入力されるステップと、
２）前記投影画像に観測されるフォトンに関する物理過程が確率分布で定義されるステップと、
３）被計測対象物の動きに関する第１の確率分布が、被計測対象物が時間と共に滑らかに変化するとし、その時間的に連続する度合いが定義されるパラメータによって表現される確率分布で定義されるステップと、
４）観測に関する第２の確率分布が、再構成画像からどのような投影像が観測されると期待されるかについて記述した確率分布で定義されるステップと、
５）上記の第１の確率分布と第２の確率分布の両方に依存して、事後確率の期待値推定（ベイズ推定）もしくは事後確率の最大化推定（ＭＡＰ推定）により画像再構成及び組織の推定が行われるステップと、
を備えたことを特徴とする動き追跡Ｘ線ＣＴ画像処理方法。
【請求項７】
被計測対象物のＸ線吸収係数に関する事前知識を設定し統計推定を行うＸ線ＣＴ画像処理装置であって、
被計測対象物が時間と共に滑らかに変化するとし、前記事前知識を第１の確率分布で設定する手段と、
前記被計測対象物のＸ線吸収係数に対してどのようなＸ線投影像が観測されるのが尤もらしいかを第２の確率分布で設定する手段と、
上記の第１の確率分布と第２の確率分布の両方に依存した統計推定を行う統計推定手段と、
を備えることを特徴とする動き追従Ｘ線ＣＴ画像処理装置。
【請求項８】
１）前記第１の確率分布が、
再構成画像の画素毎に定義されるパラメータであって、
再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる確率分布で表現され、
２）前記第２の確率分布が、
再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した確率分布で表現される、
ことを特徴とする請求項７に記載の動き追従Ｘ線ＣＴ画像処理装置。
【請求項９】
前記第１の確率分布は、再構成画像の隣接画素の空間的相関によって更に表現されたことを特徴とする請求項７に記載の動き追従Ｘ線ＣＴ画像処理装置。
【請求項１０】
前記第１の確率分布は、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することによって更に表現されたことを特徴とする請求項７に記載の動き追従Ｘ線ＣＴ画像処理装置。
【請求項１１】
請求項６の動き追従Ｘ線ＣＴ画像処理方法の各ステップを、
コンピュータに実行させるための動き追従Ｘ線ＣＴ画像処理プログラム。
【請求項１２】
請求項１１の動き追従Ｘ線ＣＴ画像処理プログラムが搭載されたＸ線ＣＴ画像処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、被計測対象物の動きに追従できるＸ線ＣＴアルゴリズム、そのアルゴリズムを有するプログラムおよび該プログラムが搭載されたＸ線ＣＴ装置に関する。
【背景技術】
【０００２】
Ｘ線ＣＴは、非侵襲・非破壊で内部構造を知ることができる強力な手法として医療用途、欠陥検査などの産業用途など幅広く用いられている。このＸ線ＣＴによる内部構造の再構成(ＣＴ画像)は、通常、内部構造を知りたい対象物が静止していることを仮定して行われてきた。しかしながら、医療用のＸ線ＣＴで対象となるもののうち、心臓・血管や肺、乳幼児などに関しては動きを完全に止めることは困難である。このような動きのある物体に対して、従来通り、静止していることを想定したＸ線ＣＴアルゴリズムを適用すると、モーションアーティファクトとよばれる偽像が生じることが知られている。
日本では、ここ１０年ほど、心筋梗塞、狭心症などの心疾患が死因の２位を占めており、アメリカでは心疾患が死因の１位であるなど心疾患は国内外で対応が求められている疾患であり、その最も有効な対応策の一つが心疾患の異常の早期発見にあるとされるため、特に心臓のＸ線ＣＴにおける正確なＣＴ画像再構成が求められている。【０００３】
従来から、モーションアーティファクトを低減させるための技術は、その応用上の重要性から主に心臓を対象として提案されている。例えば、図１４に示すように、心臓を対象としたＸ線ＣＴにおいて、撮像対象を高速にスキャンすることによって撮像時間を短くし、モーションアーティファクトを低減する先行技術１（例えば非特許文献１を参照）が知られている。また、図１５に示すように、心電図に基づく心位相の同期によりＣＴ画像を再構築する先行技術２（例えば非特許文献２〜４を参照）が知られている。さらに、心臓の動き推定に基づく動き修正することによりＣＴ画像を再構築する先行技術３（例えば非特許文献５〜９を参照）が知られている。【０００４】
一方、Ｘ線は強力であれば、再構成画像のＳＮ比が高いことから、Ｘ線被曝量の低減とＳＮ比はトレードオフの関係がある。より低いＸ線被曝量で従来と同程度のＳＮ比の再構成画像を得る技術とは、従来と同じＸ線強度，撮像枚数でも解像度の高い画像が得られるという解像度向上を図る技術や、従来と同じＸ線強度や解像度でも、より撮像枚数を減らして撮像時間を短縮できる撮像時間の短縮を図る技術や、従来と同じ解像度でもＸ線強度を弱めて被曝量を減らせる被曝量の低減を図る技術、従来と同じＸ線強度や解像度でも再構成画像に含まれるノイズを低減できるノイズ除去を図る技術、従来と同じＸ線強度や解像度でも再構成画像のコントラスト分解能（再構成画像を何階調で表せるか、すなわち、Ｘ線吸収係数の精度）を向上するコントラスト分解能の向上を図る技術がある。【０００５】
また、上記の従来技術では、統計推定を行うものと統計推定を行わないものに大きく２つに分類できる。統計推定を行わない従来技術としては、画像再構成法の主流であるフィルター補正逆投影法（ＦＢＰ：Filtered Back Projection）が知られている。ＦＢＰは、撮像対象に対して色んな方向からＸ線を照射して、その投影像(Ｘ線の吸収像、投影像のセットをシノグラムと呼ぶ)を得た後、得られた投影像を逆方向に投射して画像再構成する際に、シノグラムに周波数領域上でフィルターをかけボケ除去を行う画像再構成法である。【０００６】
一方、ＣＴ画像のように撮像対象の物体を異なる位置・角度から撮影した投影画像を複数枚用意できる場合、それら複数の投影画像の情報を処理することで、元の物体の断面像を再構成できることが知られている。投影画像を使用して元の物体の断面像を復元する画像再構成法のうち統計推定を行なう従来技術としては、ＭＡＰ（maximum a posteriori）推定法，ベイズ推定法が知られている。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】H. I. Goldberg etal., “Evaluation of ultrafast CT scanning of the adult abdomen,” Invest. Radiol., 24,537-543, 1989.
【非特許文献２】C. C. Morehouse et al., ‘‘Gated cardiac computed tomography with amotion phantom,’’ Radiology 134, 134-137, 1980.
【非特許文献３】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.
【非特許文献４】B. Ohnesorge et al., ‘‘Cardiac imaging by means ofelectrocardiographically gated multisection spiral CT: initial experience,’’ Radiology217, 564-571 , 2000.
【非特許文献５】C. J. Ritchie et al.,‘‘Correction of computed tomography motionartifacts using pixel specific back-projection,’’ IEEE Trans. Med. Imaging 15,333-342, 1996.
【非特許文献６】G. Wang et al., “A knowledge-based cone-beam x-ray CT algorithm fordynamic volumetric cardiac imaging,” Med Phys. 29(8), 1807-1822, 2002.
【非特許文献７】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.
【非特許文献８】AA. Isola et al., “Fully automatic nonrigid registration-based localmotion estimation for motion-corrected iterative cardiac CT reconstruction,” MedPhys. 37(3), 1093-1109 2010.
【非特許文献９】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.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
上記の先行技術１のように、Ｘ線源を高速に回転させる場合、遠心力が大きくなり物理的な力に耐えるように設計することが困難であり、その回転速度の高速化は頭打ちとなってきているのが現状である。現在は半周で約０．２５秒(１周約０．５秒)の速度が達成できているものの、心拍は約１秒に１周期であり、約０．２５秒の時間に心臓の体積は大きく変わることから、先行技術１のように、心臓部分が動かないと想定するにはどうしても無理が生じるという問題がある。【０００９】
また、上記の先行技術２のように、心臓のＸ線ＣＴを行う際に心電図を同時に撮り、心拍が同位相となるタイミングの投影像のみを用いて特定の位相のみの画像再構成を行う技術の場合、特定のタイミングにおける投影像のみを用いるために、結果的に被曝線量が通常よりも多くなってしまうといった問題がある。【００１０】
また、上記の先行技術３のように、心電の撮れるタイミングで同一の心位相の投影像を取得し、取得投影像から画像再構成したＸ線ＣＴ画像をもとに、他の心位相でのＣＴ画像を動きがあるものとして推定する技術の場合、心臓の一周期の長さは変わることを許容するものの、一周期内では毎回、同じ体積変化を経るものと仮定して任意の時刻における心位相を特定する。しかしながら、先行技術３では、実際には動きの早い部分などでは正確に位相を合わせることが困難であり、そもそも心臓が毎回、同じようには動いていないなどの問題がある。また、血管等は心拍ごとに位置ズレが生じる場合もあり、特に心拍が短いときと長い時とでは心臓の体積変化は異なる経過を経ることが知られており、先行技術３では、これらの要因によるアーティファクトが生じるという問題がある。また、呼吸を止めることのできない患者に対しては呼吸による動きが生じてしまうためモーションアーティファクトを防げないという問題も存在する。【００１１】
そこで、本発明は、先行技術を用いた画像再構成と比べて、動きとＣＴ画像を同時に推定して、モーションアーティファクトを低減できる、動き追従Ｘ線ＣＴ画像処理方法および動き追従Ｘ線ＣＴ画像処理装置を提供することを目的とする。
【課題を解決するための手段】
【００１２】
上記状況に鑑みて、本発明の動き追従Ｘ線ＣＴ画像処理方法は、被計測対象物の動きとＸ線吸収係数に関する事前知識をおいて統計推定を行うＸ線ＣＴ画像処理方法であって、被計測対象物が時間と共に滑らかに変化することを仮定し、動きに関する第１の確率モデル(全時刻の被計測対象物の事前知識)と観測に関する第２の確率モデル(全時刻の投影像の観測モデル)を定義するステップと、上記の第１の確率モデルと第２の確率モデルの両方に依存した統計推定を行うステップと、を備える。【００１３】
かかる構成によれば、最初から統計推定を行う確率モデルとして、被計測対象物の動きに関する第１の確率モデル(全時刻の被計測対象物の事前知識)と観測に関する第２の確率モデル(全時刻の投影像の観測モデル)を定義して、それらの確率モデルに基づいた推論を行なうことで、心電図などの同期を必要とせずに動きとＣＴ画像を同時に推定することができる。
統計推定を行う確率モデルに基づくことにより、統計の分野で確立された推論手法が適用でき、またアルゴリズムの収束性が保証されるなど、優れた最適化が行いやすい。
すなわち、従来の技術においては、得られたＣＴ再構成画像のみから動きに関するモデル（従来法は決定論的モデル）に依存して動きを推定しており、また得られた動きから観測モデルのみに依存してＣＴ再構成画像を推定しているのに対して、本発明では、動きに関する確率的モデルおよび観測に関する確率モデルの両方を考慮して、ＣＴ再構成画像と動きの両方を同時に推定することでより精度の高い推定が行える。
また、統計推定を行う確率モデルを用いているために、他の事前知識も容易に組み込むことができ、画像再構築の品質を向上できる。【００１４】
ここで、第１の確率モデル(全時刻の被計測対象物の事前知識)は、再構成画像の画素毎に定義されるパラメータであって、再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第１の確率分布で表現され、また、第２の確率モデル(全時刻の投影像の観測モデル)は、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第２の確率分布で表現されることが好ましい。【００１５】
第１の確率モデルにおいて、再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第１の確率分布で表現されるとは、異なる時刻で異なる再構成画像が仮定され、それら異なる時刻の再構成画像が時間的に滑らかに変形して得られることが確率分布の形で表現されることである。
ここで、好ましくは、再構成画像の隣接画素の空間的相関を表現する項を加えることによって、上記の第１の確率分布をさらに特徴付ける。
また、好ましくは、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第１の確率分布をさらに特徴付ける。【００１６】
また、第２の確率モデルにおいて、各時刻の再構成画像が仮に既知であると仮定した時にその再構成画像からどのような投影像が観測されると期待されるかについての情報が確率分布の形で表現される。【００１７】
また、本発明の動き追従Ｘ線ＣＴ画像処理方法において、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第１の確率分布をさらに特徴付けることが好ましい。
さらに好ましくは、本発明の動き追従Ｘ線ＣＴ画像処理方法において、再構成画像の画素毎に与えられる隠れ状態変数によって、組織クラスを表現し、上記の第１の確率分布をさらに特徴付ける。すなわち、再構成画像をさらに特徴付けるために組織クラスを表す確率変数によって、第１の確率分布をさらに特徴付ける。組織クラスを表す確率変数は、再構成画像の画素毎に、例えば、人体などの撮像対象に関し、どのような組織（筋肉などの通常細胞、脂肪などの軟細胞、骨など）がどのようなＸ線吸収係数を有するかを表現する。再構成画像は時間的に滑らかに変形することが第１の確率分布に表現されるが、その表現は各時刻の再構成画像が時間的に直接相関をもつことを仮定する確率分布による表現の他、各時刻の組織クラスが時間的に滑らかに変位すること仮定することで組織クラスを通じた間接的な表現を行うことができる。また、これらの組織がどの程度の割合で分布するか、また、それぞれの組織がどの程度、空間的に連続して分布しやすいか、の分布情報が確率分布の形で表現される。
なお、組織分布に関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と組織推定の精度のさらなる向上を図ることができる。【００１８】
また、本発明の動き追従Ｘ線ＣＴ画像処理方法において、画素毎に所属する組織クラスを表現する隠れ状態変数によって、第１の確率分布をさらに特徴付けることが好ましい。
画素毎に所属する組織クラスの分布を制御する隠れ状態変数により、事前に想定した組織クラスのＸ線吸収係数の想定値が実際の組織クラスのＸ線吸収係数から大きくズレたり、想定していなかった組織クラスが存在していた場合に影響を受け難くすることができる。また、この隠れ状態変数を統計推定するにより、組織依存のＣＴ値が存在するという事前知識がどの程度、推定される解である再構成画像へ影響するか、という影響度の度合いを制御することができる。【００１９】
ここで、本発明の動き追従Ｘ線ＣＴ画像処理方法において、上記の組織クラスの事前分布としてボルツマン分布を用いることが好ましい。
各組織の事前分布としてボルツマン分布を用いることとした理由は、被計測対象物が時間と共に滑らかに変化することを仮定したことから、隣接するフレーム間でも組織は時間と共に滑らかに変化することを仮定するべきであるが、その仮定を隣接するフレーム間の組織クラスの時間相関によっての表現しやすいことに注目したものである。
その他、組織が満たすべき性質である、組織の空間的連続性や、どの組織がどの程度の割合で存在しやすいかも同時に表現できる特性をもつ。 時間相関や空間相関の強さ、標準的な存在割合を表すパラメータは、組織クラスごとに専門家による人手による調整や実際の投影データへの当てはまりの良さを基準とした学習による調整が可能である。
なお、上記の各組織のとる典型的なＸ線吸収係数の値を表現するパラメータ分布としてはガウス分布を好適に用いることができ、また、各組織の空間的に連続する度合いを表現するパラメータ分布としてはボルツマン分布を好適に用いることができる。組織クラスの分布を制御する隠れ状態変数を導入した際には、各組織のＸ線吸収係数を表現するパラメータ分布はガウス分布の重み付き和である混合ガウス分布を好適に用いることができる。【００２０】
また、上記の本発明の動き追従Ｘ線ＣＴ画像処理方法における統計推定は、具体的には、事後確率を最大化する値による推定(ＭＡＰ推定)、或いは、事後確率の期待値によるベイズ推定で行う。
ここで、ＭＡＰ推定は、再構成画像の各画素の領域が属する組織とＸ線吸収係数の両方に関する事後分布において最も尤もらしい組合せを推定する。再構成画像の各画素の領域が属する組織との組合せを推定することは、組織クラス依存性を含めることになる。この組織クラス依存性を含めることで、例えば、空気の組織クラスはピクセル同士の連結がしやすい（空気中に他の組織は入らない）などの組織クラスに依存した知識を含めることができるのである。これにより、より現実の状況に沿った分布表現が可能となり結果的に解の推定精度を高めることが期待される。ＭＡＰ推定によれば、再構成画像の事後確率を求めるのに必要な積分計算が省略することができるため、最適化を速めることができ、より高速で再構成画像を求めることができる。【００２１】
また、(近似)ベイズ推定は、Ｘ線吸収係数と再構成画像の各画素の領域が属する組織の両方に関する事後分布を推定するにあたり、事後分布をよく近似可能な試験分布を用いて近似する。ベイズ推定の場合も、ＭＡＰ推定の場合と同様に、組織クラス依存性を含めることができる。これにより、より柔軟性をもった分布表現が可能となる。ベイズ推定によれば、推定の不確実性を考慮に入れた推論が可能で、またパラメータの学習が容易であり、より正確にＸ線ＣＴ画像処理を行うことができる。【００２２】
また、本発明の他の観点によれば、本発明の動き追従Ｘ線ＣＴ画像処理方法は、下記１）〜５）のステップを備える。
１）投影画像と少なくともＸ線強度を含む計測条件が入力されるステップ
２）投影画像に観測されるフォトンに関する物理過程が確率分布で定義されるステップ
３）被計測対象物の動きに関する第１の確率モデル(全時刻の被計測対象物 の事前知識)が、被計測対象物が時間と共に滑らかに変化することを仮定して、各組織が時間的に連続する度合いが定義されるパラメータによって特徴付けられる第１の確率分布の形で定義されるステップ
４）観測に関する第２の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについての情報が第２の確率分布の形で定義されるステップ
５）上記の第１の確率モデルと第２の確率モデルの両方に依存して、事後確率の期待値推定（ベイズ推定）もしくは事後確率の最大化推定（ＭＡＰ推定）により画像再構成及び組織クラス推定が行われるステップ【００２３】
上記２）のフォトンに関する物理過程が確率分布は、具体的には、投影画像に観測されるフォトンに関する確率モデルがポアソン分布などで表現される。これにより、より現実の観測過程に近い物理過程を表現することができ、フォトンノイズなどによる観測の不確定性を表現できる。
上記３）〜５）についての説明は、上述の説明と同様である。【００２４】
次に、本発明の動き追従Ｘ線ＣＴ画像処理装置について説明する。
本発明の動き追従Ｘ線ＣＴ画像処理装置は、被計測対象物の動きとＸ線吸収係数に関する事前知識をおいて統計推定を行うＸ線ＣＴ画像処理装置であって、被計測対象物が時間と共に滑らかに変化することを仮定し、動きに関する第１の確率モデル(全時刻の被計測対象物の事前知識)と観測に関する第２の確率モデル(全時刻の投影像の観測モデル)を定義するモデル定義手段と、上記の第１の確率モデルと第２の確率モデルの両方に依存した統計推定を行う統計推定手段を備える。【００２５】
ここで、第１の確率モデル(全時刻の被計測対象物の事前知識)は、再構成画像の画素毎に定義されるパラメータであって、再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第１の確率分布で表現され、また、第２の確率モデル(全時刻の投影像の観測モデル)は、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第２の確率分布で表現されることが好ましい。【００２６】
第１の確率モデルにおいて、再構成画像が時間的に連続する度合いを表現するパラメータによって特徴付けられる第１の確率分布で表現されるとは、異なる時刻で異なる再構成画像が仮定され、それら異なる時刻の再構成画像が時間的に滑らかに変形して得られることが確率分布の形で表現されることである。
ここで、好ましくは、再構成画像の隣接画素の空間的相関を表現する項を加えることによって、上記の第１の確率分布をさらに特徴付ける。
また、好ましくは、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第１の確率分布をさらに特徴付ける。【００２７】
また、第２の確率モデルにおいて、各時刻の再構成画像が仮に既知であると仮定した時にその再構成画像からどのような投影像が観測されると期待されるかについての情報が確率分布の形で表現される。【００２８】
また、本発明の動き追従Ｘ線ＣＴ画像処理装置において、上述の本発明の動き追従Ｘ線ＣＴ画像処理方法と同様、再構成画像の各画素の値が組織に依存して特定の値をとりやすいことを表現する確率分布と組織も同様に時間的に連続に変位することを表す確率分布によって、上記の第１の確率分布をさらに特徴付けることが好ましい。
さらに好ましくは、本発明の動き追従Ｘ線ＣＴ画像処理装置において、再構成画像の画素毎に与えられる隠れ状態変数によって、組織クラスを表現し、上記の第１の確率分布をさらに特徴付ける。すなわち、再構成画像をさらに特徴付けるために組織クラスを表す確率変数によって、第１の確率分布をさらに特徴付ける。組織クラスを表す確率変数は、再構成画像の画素毎に、例えば、人体などの撮像対象に関し、どのような組織（筋肉などの通常細胞、脂肪などの軟細胞、骨など）がどのようなＸ線吸収係数を有するかを表現する。再構成画像は時間的に滑らかに変形することが第１の確率分布に表現されるが、その表現は各時刻の再構成画像が時間的に直接相関をもつことを仮定する確率分布による表現の他、各時刻の組織クラスが時間的に滑らかに変位すること仮定することで組織クラスを通じた間接的な表現を行うことができる。また、これらの組織がどの程度の割合で分布するか、また、それぞれの組織がどの程度、空間的に連続して分布しやすいか、の分布情報が確率分布の形で表現される。
なお、組織分布に関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と組織推定の精度のさらなる向上を図ることができる。【００２９】
また、本発明の他の観点によれば、本発明の動き追従Ｘ線ＣＴ画像処理装置は、下記１）〜５）の手段を備える。
１）投影画像と少なくともＸ線強度を含む計測条件を入力する手段
２）投影画像に観測されるフォトンに関する物理過程を確率分布で定義する手段
３）被計測対象物の動きに関する第１の確率モデル(全時刻の被計測対象物 の事前知識)が、被計測対象物が時間と共に滑らかに変化することを仮定して、その時間的に連続する度合いを定義するパラメータによって特徴付けられる第１の確率分布を定義する手段
４）観測に関する第２の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時に、その再構成画像からどのような投影像が観測されると期待されるかについて記述した第２の確率分布を定義する手段
５）上記の第１の確率モデルと第２の確率モデルの両方に依存して、事後確率の期待値推定（ベイズ推定）もしくは事後確率の最大化推定（ＭＡＰ推定）により画像再構成及び組織クラス推定を行う手段【００３０】
上記２）のフォトンに関する物理過程が確率分布は、具体的には、投影画像に観測されるフォトンに関する確率モデルがポアソン分布などで表現される。これにより、より現実の観測過程に近い物理過程を表現することができ、フォトンノイズなどによる観測の不確定性を表現できる。【００３１】
また、本発明の動き追従Ｘ線ＣＴ画像処理プログラムは、被計測対象物の動きとＸ線吸収係数に関する事前知識をおいて統計推定を行うＸ線ＣＴ画像処理プログラムであって、コンピュータに、下記１）〜５）のステップを実行させる。
１）投影画像と少なくともＸ線強度を含む計測条件を入力するステップ
２）前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップ
３）被計測対象物の動きに関する第１の確率モデル(全時刻の被計測対象物の事前知識)が、各時刻の再構成画像が時間的に連続する度合いが定義されるパラメータによって特徴付けられる第１の確率分布で定義するステップ
４）観測に関する第２の確率モデル(全時刻の投影像の観測モデル)が、再構成画像が仮に既知であると仮定した時にその再構成画像から、どのような投影像が観測されると期待されるかについて記述した第２の確率分布で定義するステップ
５）上記の第１の確率モデルと第２の確率モデルの両方に依存して、事後確率の期待値推定（ベイズ推定）もしくは事後確率の最大化推定（ＭＡＰ推定）により画像再構成及び組織クラス推定を行うステップ【００３２】
また、本発明は、上記のＸ線ＣＴ画像処理プログラムを搭載したＸ線ＣＴ画像処理装置を提供することができる。
【発明の効果】
【００３３】
本発明の動き追従Ｘ線ＣＴ画像処理方法および動き追従Ｘ線ＣＴ画像処理装置によれば、動きとＣＴ画像を同時に推定して、従来技術を用いた画像再構成と比べて、モーションアーティファクトを低減できるといった効果を有する。
また、被計測対象物が時間的に滑らかに変形するという仮定に基づいた推定方法であるため、心電図などによる同期を取る必要がなく、心電図などの同期が取れない予測不可能な体動や内臓運動に対応が可能な他、全時刻の投影像を利用した推定が行えるため被曝量の低減が期待される。
【図面の簡単な説明】
【００３４】
【図１】通常のＣＴモデルと動きを想定したＣＴモデルの説明図
【図２】隣接フレーム間における画素の動きの説明図
【図３】本発明の動き追従Ｘ線ＣＴ画像処理方法のフロー図（１）
【図４】本発明の動き追従Ｘ線ＣＴ画像処理方法のフロー図（２）
【図５】組織クラスごとのＸ線吸収係数の分布図
【図６】動きと組織クラス、組織クラスの事前分布を制御する隠れ状態変数を想定したＣＴの確率モデルにおける確率変数間の依存関係を表す説明図
【図７】隠れ状態変数Ｂの導入効果の説明図
【図８】時間と共に変形する物体の説明図
【図９】時間と共に変形する物体に対するＣＴ画像再構成の結果を示す図
【図１０】心臓の動きを伴う胸部モデルの説明図
【図１１】心臓の動きを伴う胸部モデルに対するＣＴ画像再構成の結果を示す図
【図１２】組織クラスの分布を制御する隠れ状態変数ｂの導入効果を示す図
【図１３】ＣＴ画像再構成の結果を示す図
【図１４】先行技術１の説明図
【図１５】先行技術２の説明図
【発明を実施するための最良の形態】
【００３５】
以下、本発明の実施形態について、図面を参照しながら詳細に説明していく。ただし、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更または変形が可能である。【００３６】
本発明の動き追従Ｘ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、動きとＸ線吸収係数に関する事前知識をおいて統計推定を用いて画像の再構成を行う。
図１を参照し、通常のＣＴモデルと動きを想定したＣＴモデルの説明を行う。図１（１）に示すように、従来のフィルター補正逆投影法（ＦＢＰ）のＣＴモデルの場合、動きのある心臓などの被計測対象物であっても、対象は静止していると仮定して１つの対象(ＣＴ像)から投影が得られると考えていた。本発明の場合、図１(２)に示すように、被計測対象物は初めから動いているものと仮定し、投影はこれらの異なる対象から得られるものと考える。【００３７】
そこで、静止した被計測対象物、あるいは周期的に変動する被計測対象物に対して投影を行っているという仮定はおかず、被計測対象物は時間と共に滑らかに変化することを仮定し、この仮定をもとに確率モデルを構築する。
被計測対象物は時間と共に滑らかに変化することを仮定するとは、ＣＴ画像の画素毎に動き（変形）パラメータを設け、画素毎に隣接フレーム間での動き（変形）を推定することである。具体的には、図２に示すように、ｊ番目の画素の異なる時刻間の再構成画像（ｋ−１番目のフレームとｋ番目のフレーム）の間での、動きをベクトルｍｊ（ｋ）で表す。
以下では、時間と共に滑らかに変化する被計測対象物に対して投影を行っていることを統計モデルによって記述し、この統計モデルに基づきベイズ推定を行うことで画像再構成できることを説明する。【００３８】
図３と図４に、本発明の動き追従Ｘ線ＣＴ画像処理方法の実施形態のフローを示す。
本発明の動き追従Ｘ線ＣＴ画像処理方法のフローの一実施形態としては、図３のフローに示すように、投影画像とＸ線強度などの計測条件を入力し（ステップＳ１１）、投影画像の観測フォトンに関して物理モデルで表現し（ステップＳ１２）、Ｘ線吸収係数に関する事前知識を、被計測対象物が時間と共に滑らかに変化することを仮定し、再構成画像が時間的に連続する度合いを表現するパラメータで特徴付けた第１の確率モデル(全時刻の被計測対象物の事前知識)と第２の確率モデル(全時刻の投影像の観測モデル)とで表現し（ステップＳ１３）、事後確率の期待値推定（ベイズ推定）で、動き，組織クラスを推定および面像再構成する（ステップＳ１４）。【００３９】
また、本発明の動き追従Ｘ線ＣＴ画像処理方法のフローの他の実施形態としては、図４のフローに示すように、画素毎に所属する組織クラスを表現する隠れ状態変数によって、上記の第１の確率分布をさらに特徴付ける（ステップＳ２４）。
画像再構成は不良設定性を有しているため、適当な事前知識を用いて解である再構成画像(X線吸収係数)に何らかの制約を課すことで解決する。この制約のうち、本発明の肝となるのが隣接するフレーム間の再構成画像に与えられる時間的制約である。非計測対象物が時間的に滑らかに変形することを想定しているため、図２のとおり隣接するフレーム間の再構成画像も各画素jごとに動きｍｊ（ｋ）に応じた動きを生じる。したがって、ｋ−１番目のフレームにおける再構成画像を動きｍｊ（ｋ）にしたがって変形させたものは、ｋ番目のフレームにおける再構成画像と強い時間相関をもつべきである。
この再構成画像に対する時間的制約は、再構成画像間の直接の時間相関と、後述の再構成画像の背後に潜在的に存在すると仮定される物質・組織を介した間接的な時間相関の両方、もしくは一方によって確率分布の形で表現される。【００４０】
物質・組織を介した制約は以下のとおりでる。図５に示すように、人体は脂肪，水、血、筋肉，骨，金属（例えばインプラント）といった限られた物質・組織から構成され、それぞれの物質・組織(以下、単に組織とよぶ)がある典型的なＸ線吸収係数をとりやすいといった制約を確率分布の形で特徴づけることができる。また、被計測対象の動きに沿って、組織が変形することを隣接するフレーム間の再構成画像の時間相関と同様な形で表現することができ、この制約が確率分布の形で特徴付けられる。その他、それぞれの組織は空間的な連続性を有しており、各組織が人体に占める割合も大まかに分かっているといった知識も組織に対する制約として同様に確率分布の形で特徴付けられる。【００４１】
図６は、図１(２)を本発明の確率モデル通りに各変数間の因果関係を具体的に書いたものである。スキャン時間を十分短い時間幅のフレームに分け、各フレーム内では対象は同一のＣＴ画像で表されるとする（図６では、これを、ｘ（１），ｘ（２），・・・のように表す。)。 実際のスキャンでは、０．３秒ほどの短い時間に数百から数千の投影が得られるので、各フレームのｘからは複数の投影ｙが得られる。不良設定性に対処するため、組織クラスＺを導入し、Ｚが与えられたもとで、Ｘ線吸収係数ｘはある程度決まった値を取るとする。また、Ｂは事前分布の形を制御する隠れ状態変数で、事前分布をよりロバストなものにしている。【００４２】
図７は、隠れ状態変数Ｂを導入することによる効果を示す図である。事前分布の形を制御する隠れ状態変数Ｂがない場合、例えば事前分布の分散を小さくして曖昧性を低くした場合は、それに依存して事前知識の与える効果も強くなってしまう。Ｂを導入することにより事前分布を混合分布で表すことができ、事前知識の曖昧さと効果の強さを別々に制御することができる。【００４３】
そして、事後分布に基づいた統計推定を行う。再構成画像に関する事後分布が最大となる値(ＭＡＰ推定)もしくは期待値(ベイズ推定)を求めることによって画像再構成が行える。一般に事後確率の期待値を求めるベイズ推定のほうが擾乱に強い推定が行えるが、事後確率の計算には高次元の隠れ状態変数に関する積分計算が含まれるため、解析的に実行することは困難であることから、事後確率を近似的に求める近似手法を適用することでこの計算困難さを克服する。なお、事後分布推定による統計推定においては、隠れ状態変数に関する積分計算を行わず隠れ状態変数と再構成画像の両方を含んだ同時事後確率の最大化によるＭＡＰ推定を適用することでこの計算困難さを克服しても構わない。【００４４】
（問題の定式化）
先ず、Ｘ線ＣＴ画像のように、様々な方向から投影されたＴ個の投影データを、Ｄ＝｛Ｙ（１），・・・，Ｙ（Ｔ）｝と表すことにする。各々のデータＹ（ｔ）は、ｔ番目の投影によって検出器で検出されるデータの集合であり、Ｙ（ｔ）＝｛ｙ１（１），・・・，ｙＩ（ｔ）｝となる。ｙｉ（ｔ）は、ｉ番目の検出器で検出された光子の数を表す。
短時間に十分、多くの投影が得られているときには、その短時間内では、Ｘ線ＣＴ画像が同一であると想定したほうが推定すべきＸ線ＣＴ画像を減らすことができ計算量を抑えることができる。そこで、各投影をＫ通りのフレームに分け、各フレーム内では同一のＸ線ＣＴ画像から投影像が得られると仮定する。ｋ番目のフレームに所属する投影データのインデックスの集合をＳ（ｋ）と表し、ｋ番目のフレームにおけるＸ線ＣＴ画像をＸ（ｋ）と表すことにする。Ｘ線は物質を透過する際に指数的に減衰することから、観測されると期待されるＸ線光子の平均ｙｉ（ｔ）は、下記数式１のように表すことができる。【００４５】
【数１】
0006044046.tif 000002
【００４６】
ここで、ｘｊ（ｋ）はｋ番目のフレームにおける撮像対象のＸ線吸収係数をラスタスキャンして得られるＪ次元ベクトルｘ＝｛ｘ１，・・・，ｘＪ｝のｊ番目のピクセルのＸ線吸収係数である。ｖは、Ｘ線源から放出される光子数（物体が何も置かれていないときに観測され得る光子数）を表す。また、ｌｉｊ（ｔ）は角度θ（ｔ）から投影された際のｉ番目の検出器で検出される投影線とｊ番目のピクセルが交差する距離に比例するものであり、ｌｉｊ（ｔ）xｊがｔ番目の投影のｉ番目の検出器に入射するＸ線に対するｊ番目のピクセル領域の実効的なＸ線吸収係数を表す。【００４７】
次に、本発明の動き追従Ｘ線ＣＴ画像処理方法および動き追従Ｘ線ＣＴ画像処理プログラムにおいて、投影画像に観測されるフォトンに関する物理モデル、事前知識としての組織に関する事前分布、すなわち、撮像対象の各組織の存在割合を表現するパラメータ分布、各組織がとりやすい典型的なＸ線吸収係数を表現するパラメータ分布、各組織の空間的に連続する度合いを表現するパラメータ分布、各組織の時間的に連続する度合いを表現するパラメータ分布について、以下説明を行う。【００４８】
（観測されるフォトンに関する物理モデル）
Ｘ線ＣＴにおける観測データの主なノイズの要因は、検出される光子（フォトン）のショットノイズであると考えられる。そこで、光子（フォトン）の観測データは、投影毎、検出器毎に独立なポアソン分布に従って生成されるとしてモデル化する。
なお、ポアソン分布に特に限定されるものではなく、物理過程をより良く表現できる他の物理モデルも適用可能である。
上述の数式１は、ポアソン分布の平均を表しており、数式１の結果を用いて、下記数式２で表わす。ノイズは投影、画素（ピクセル）毎に独立であることから、数式２は全体として独立なポアソン分布の積で表されている。数式２において、Ｓ（ｋ）は、ｋ番目のフレームに含まれる投影のインデックスの集合を表している。【００４９】
【数２】
0006044046.tif 000003
【００５０】
（動きに関する事前分布）
動きθが与えられたもとでの組織クラスの分布はボルツマン分布で表されると仮定して、組織クラスに関する事前分布を下記数式３に示すようにモデル化する。ここで、Ａ１は正規化定数であり、Ｈ１（Ｚ）はどのような組織クラスの組み合わせが尤もらしいか表すエネルギー関数である。エネルギー関数Ｈ１（Ｚ）は、下記数式４で表される。
ここで、ｚｊ(ｋ) = [ｚｊ１(ｋ), ・・・, ｚｊＣ(ｋ) ] ∈ {０，１}Ｃは、ｋ番目のフレームにおいてｊ番目のピクセルが属する組織クラスを表すベクトルであり、組織クラスｃに属する場合はｚjc(ｋ) ＝１、その他の成分は０となる。大文字のＺは小文字の変数ｚｊ(ｋ)の集合を表す。【００５１】
【数３】
0006044046.tif 000004
【００５２】
【数４】
0006044046.tif 000005
【００５３】
上記数式４において、η（ｊ）は、ｊ番目のピクセルの４近傍、Ｎ（ｊ）はｊ番目のピクセルを含むその周囲９×９（ピクセル）。１項目は組織の割合、２項目は同一フレームにおける組織の空間的なめらかさ(つながりやすさ)、３項目は隣接するフレーム間での組織の時間的になめらかな変位・変形を表現するものである。また、数式４中、Ｇ（ｉ，ｊ，ｍ（ｋ）（θ））はＨ１（Ｚ）の３項目の重みであり、下記数式５で表されるもので、隣り合うフレーム(ｋ−１番目とｋ番目) において、ｋ−１番目のフレームの各ピクセルが現在の動きｍ（ｋ）（θ）にしたがって動いたとして、ｋ番目のフレームの各ピクセルとの距離を計算し、その距離が近いほど重みが大きくなる。【００５４】
【数５】
0006044046.tif 000006
【００５５】
（再構成画像Ｘに関する事前分布）
再構成画像Ｘの事前分布を下記数６に示す。組織クラスＺと隠れ状態変数Ｂが与えられたもとで、Ｘ線吸収係数はガウス分布として表現している。以下に説明するように、ガウス分布の平均は組織クラスごとに設定する。また、その組織クラスで典型的なＸ線吸収係数の値への接近の度合いは隠れ状態変数ｂによって制御される。各フレームｋ、各画素ｊにおいてｂj１（ｋ）, ｂj２（ｋ）はいずれか一方が１をとりいずれか一方が０をとる変数である。組織クラスｃごとに与えられるその組織クラスｃに典型的なＸ線吸収係数の値ｖｃ１, ｖｃ２は通常、同じ値ｖｃ１＝ｖｃ２＝ｖｃが設定される。たとえば、平均ｖｋの値は、各々の組織（空気，脂肪，筋肉，骨，金属）ごとにそれぞれ、ｖ１＝0.000，ｖ２＝0.018，ｖ３＝0.022，ｖ４＝0.040，ｖ５＝0.120と設定する。また、分散パラメータεｃ１２, εｃ２２は組織クラスｃごとに固定した値であるが、これは一方が大きく、一方が小さい正の値をとるパラメータである。ｖｃ１, ｖｃ２,εｃ１２, εｃ２２はいずれも専門家による人手あるいは実際に得られた多くの投影像への当てはまりの良さを基準としたパラメータ学習によって与えられる。
数式６の１項目は典型的なＸ線吸収係数の値への再構成画像が表すＸ線吸収係数の近接度の度合いを決定する。隠れ状態変数ｂj１（ｋ）, ｂj２（ｋ）はいずれか一方が１、いずれか一方が０をとる変数であるため、再構成画像の組織クラスｃの典型的なＸ線吸収係数への値の近接度は組織クラスｃのみによって決定されるのではなく各フレームｋ、各画素ｊごとにｂj１（ｋ）, ｂj２（ｋ）のとる値によっても制御されることになる。このような組織クラスｃごとに固定された分散パラメータεｃ１２, εｃ２２の他に、与えられた投影像から推定される隠れ状態変数ｂj１（ｋ）, ｂj２（ｋ）を用意し推定することで各フレーム、各画素ごとに適応的に近接度を制御可能とし、
これによって組織クラスごとにどのＸ線吸収係数をとりやすいかに関する不確かさとその影響の強さを別々に制御することができる。数式６の２項目と３項目はそれぞれ再構成画像が空間的に滑らかに変化すること、動きベクトルｍに基づいて時間的に滑らかに変形することを表す項であり、ＧｘspatialとＧｘteemporalは空間相関の強さと時間相関の強さを制御するパラメータである。【００５６】
【数６】
0006044046.tif 000007
【００５７】
ここで、隠れ状態変数Ｂはボルツマン分布に従うと仮定した場合（下記数式７）、隠れ状態変数Ｂのどのような組み合わせが尤もらしいかを表すエネルギー関数Ｈ２（Ｂ）は下記数式８で表される。ここで、数式８の１項目は２値の隠れ変数の集合であるＢがそれぞれ２値のいずれの値をとりやすいかを表し、２項目は分布の空間的なめらかさ（隣接するピクセルでの同じ分布の取りやすさ）を制御するものである。【００５８】
【数７】
0006044046.tif 000008
【００５９】
【数８】
0006044046.tif 000009
【００６０】
また、ｋ−１番目のフレームからｋ番目のフレームへのｊ番目のピクセルの動きを表すｍj（ｋ）（θ）は、下記数式９で表される。数式９は、少数の代表動きベクトルθi（ｌ）(ｉ＝１,・・・,Ｊ, ｌ＝１,・・・，Ｋ)からｋ番目のフレームにおけるｊ番目のピクセルの動きベクトルｍj（ｋ）（θ）を補間により求める式である。空間的、時間的に近い代表動きベクトルθi（ｌ）(ｉ＝１,・・・,Ｊ, ｌ＝１,・・・，Ｋ)の重みが大きくなるようになっている。なお、数式９において、ｒｉはｉ番目の画素の位置ベクトルを表す。  【００６１】
【数９】
0006044046.tif 000010
【００６２】
（事後分布のベイズ推定）
ＣＴ画像の再構成では、Ｘ線吸収係数ｘと組織クラスｚと隠れ状態変数ｂは、観測データＤと被計測対象物の動きθの事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ，θ）として推定される。事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ，θ）は、ベイズの定理により分解でき、それぞれの事前分布と尤度の積に比例する形で、下記数式１０のように表される。ＭＡＰ推定においては事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ，θ）の最大値を求めることにより、高次元の隠れ状態変数に関する和計算の必要がなくなる。【００６３】
【数１０】
0006044046.tif 000011
【００６４】
(近似)ベイズ推定では、計算困難な事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ，θ）を近似する試験分布ｑ（Ｘ，Ｚ，Ｂ）によって評価される。なお、ベイズ推定自体は、事後分布を計算しその期待値を推定値とする方法であって、事後分布やその期待値が解析的に計算可能なら試験分布は不要である。
試験分布ｑ（Ｘ，Ｚ，Ｂ）は、ＫＬ距離と呼ばれる分布の近さを表す指標を最小化するように決定される。ここでのＫＬ距離は下記数式１１で定義される。【００６５】
【数１１】
0006044046.tif 000012
【００６６】
ここで、＜・＞ｑ（Ｘ，Ｚ，Ｂ）は、試験分布ｑ（Ｘ，Ｚ，Ｂ）に関して積分計算を行うことを表す演算子である。すなわち、試験分布ｑ（Ｘ，Ｚ，Ｂ）に関して期待値を取ることを意味する。ＫＬ距離は常に非負であり、ｑ＝ｐのときにのみ０（ゼロ）になる。事後分布の計算困難性に対処するために、計算可能な分布の集合の中において事後分布に最も近いものを探索する。
試験分布ｑ（Ｘ，Ｚ，Ｂ）は、ＫＬ距離（上記数式１１）の最小化が可能であるようなもの、あるいは近似的な最小化が可能であるようなものであれば、任意に選ぶことができる。そのため、計算を簡単にするために、試験分布ｑ（Ｘ，Ｚ，Ｂ）には、下記数式１２に示すように、試験分布ｑ（Ｘ，Ｚ，Ｂ）に対して因子化仮定と呼ばれる変数ｘとｚとｂの間の独立性の仮定を行う。すなわち、Ｘ，Ｚ，Ｂは独立で、さらにフレーム、ピクセル毎にも独立であるとする。  さらに、ｘｊに関する分布ｑ（ｘｊ）がガウス分布になるものと仮定して、下記数式１３に示されるようにモデル化する。これらの仮定の下で、ｑ（Ｘ）,ｑ（Ｚ）,ｑ（Ｂ）を交互に最適化する交互最適化を用い最適試験分布を求めている。【００６７】
【数１２】
0006044046.tif 000013
【００６８】
【数１３】
0006044046.tif 000014
【００６９】
ｑ（Ｘ）とｑ（Ｂ）を現在、得られている推定値に固定し、ｑ（Ｚ）のみを上記数式１１のＤＫＬが最小化できるよう最適化する際には、解析解が得られる。同様に、ｑ（Ｘ）とｑ（Ｚ）を現在、得られている推定値に固定し、ｑ（Ｂ）のみを上記数式１１のＤＫＬが最小化できるよう最適化する際にも、解析解が得られる。
ｑ（Ｚ）とｑ（Ｂ）を現在、得られている推定値に固定しても、ｑ（Ｘ）は解析解が得られないので、上記数式１３におけるｑ（Ｘ）のパラメータμ，σは数値最適化によって求める。これらの分布ｑ（Ｘ）, ｑ（Ｚ）, ｑ（Ｂ）は、事前に決められた所定の収束条件を満たすまで最適化が繰り返され、解が求められる。【００７０】
以上、本発明の動き追従Ｘ線ＣＴ画像処理方法および動き追従Ｘ線ＣＴ画像処理プログラムにおける処理内容について説明した。以下の実施例では、変形する物体を人工的に生成し、その投影像から画面再構成を行い、従来技術による画面再構成と比較して評価したものである。下記の評価実験から、本発明の動き追従Ｘ線ＣＴ画像処理方法および動き追従Ｘ線ＣＴ画像処理プログラムの有用性の理解が深まるであろう。
【実施例１】
【００７１】
実施例１では、変形する物体を人工的に生成し、その投影像から本発明の動き追従Ｘ線ＣＴ画像処理方法により画面再構成を行う計算機シミュレーションを行った結果を説明する。
先ず、図８に示す画像は、人工的に生成した変形する物体（６４×６４ピクセル）であり、動きのある物体を模擬するものである。図８に示すように、それぞれの円形の画面の中央部分が物体を表しており、これが時刻と共に（ａ）→（ｂ）→（ｃ）→（ｄ）→（ｅ）と膨張していくものとした。物体の変形は、縦と横にそれぞれ異なる比率で膨張していくものとした。
Ｘ線ＣＴ画像の再構成の結果を図９に示す。図９は、真の画像と、動きを想定しない従来法として最もよく用いられているフィルター補正逆投影法（ＦＢＰ）による画像再構成の結果と、本発明による画像再構成（提案法）の結果とを比較したものである。
それぞれの手法による再構成画像と真の画像との差をＰＳＮＲ(Peak signal to noise ratio)によって評価した。ＰＳＮＲは下記数式１４のように定義されるものである。【００７２】
【数１４】
0006044046.tif 000015
【００７３】
上記数式１４において、ＭＳＥは真の画像と再構成画像の画素値間の平均二乗誤差を表し、ｍａｘｉｙｉは、真の画像の画素値の最大値を表す。図９の各図の下部にＰＳＮＲの値を付記している。
図９の結果から、本発明はＰＳＮＲで平均して約３．８（ｄＢ）の改善をもたらしていることがわかる。約３．８（ｄＢ）の改善は平均二乗誤差でいうと約４２％の誤差を低減していることを意味する。得られる画像からも中心の変形する物体の動きをよく表現できていることがわかる。【００７４】
次に、図１０に示す画像は、心臓の動きを伴う胸部モデルの画像であり、心臓の動きを伴う胸部を模擬するものである。図１０に示すように、心臓に対応する中心の上の楕円がだんだん大きくなるように、すなわち時刻と共に（ａ）→（ｂ）→（ｃ）→（ｄ）→（ｅ）と心臓に対応する楕円が膨張していくものとした。楕円の変形は、縦と横にそれぞれ異なる比率で膨張していくものとした。胸部モデルの画像サイズは、２５６×２５６（ピクセル）、投影数は３６０°の間に３６０枚とした。
Ｘ線ＣＴ画像の再構成の結果を図１１に示す。図１１は、真の画像と、動きを想定しない従来法として最もよく用いられているフィルター補正逆投影法（ＦＢＰ）による画像再構成の結果と、本発明による画像再構成（提案法）の結果とを比較したものである。
図１１の結果から、本発明はＰＳＮＲで平均して約１．１（ｄＢ）の改善をもたらしていることがわかる。約１．１（ｄＢ）の改善は平均二乗誤差でいうと約２２％の誤差を低減していることを意味する。得られる画像からも心臓の動きをよく表現できていることがわかる。
【実施例２】
【００７５】
本発明の動き追従Ｘ線ＣＴ画像処理方法では、上述したように、統計推定を行う確率モデルを用いているために、他の事前知識も容易に組み込むことができ、画像再構築の品質を向上できる。例えば、各組織の代表的なＸ線吸収係数を考慮に入れて推定を行うことで推定精度を向上させることができる。この際、各組織クラスの分布を制御する隠れ状態変数を含めることで推定精度の向上、ロバストな推定ができる。その効果を調べた実験について、以下説明する。【００７６】
図１２は、想定した組織クラスのＸ線吸収係数と、実際に実験で用いた物体のもつＸ線吸収係数の値を示している。ここでは、空気、脂肪、筋肉、骨の４つの組織クラスが存在すると仮定し、それぞれの組織クラスで代表的なＸ線吸収係数があるものとしている。
実際の投影では、想定した代表値からＸ線吸収係数の値のズレが生じることを考えて、図中に符号ａで示した線のように、脂肪と骨のクラスのＸ線吸収係数は１０％だけＸ線吸収係数が高いものとした。また図中に符号ｂで示した線は、その他の外れ値のとるＸ線吸収係数の値である。各組織クラスのうち、脂肪と筋肉に関して想定した代表的なＸ線吸収係数から５％大きな値をとる物体に対してＸ線ＣＴを行った。【００７７】
図１３は、上記のような条件設定で、Ｘ線ＣＴ再構成の結果を比較したものである。図１３において、（ａ）は真の画像であり、（ｂ）はＦＢＰによる再構成画像であり、（ｃ）は本発明による再構成画像で隠れ状態変数Ｂの導入が無しのものであり、（ｄ）は本発明による再構成画像で隠れ状態変数Ｂの導入が有りのものである。また（ｅ）は、推定された隠れ状態変数Ｂを表している。
隠れ状態変数Ｂを導入した場合と導入しない場合とで、ＰＳＮＲで０．８（ｄＢ）の改善（平均二乗誤差で１７％程度の改善）となった。従来法であるフィルター補正逆投影法（ＦＢＰ）と比較すると１７．６（ｄＢ）の改善となっている。また、各組織におけるＸ線吸収係数が想定からズレた場合を考えているが、５％と小さいシフトしか想定していない。このため、より大きなズレが生じる場合には、より大きな効果が得られると期待される。【００７８】
以上の結果は、比較的、確率モデルを簡単化したものでの結果であり、推定精度には向上の余地が存在するものの、推定精度自体は実用化のための十分なレベルを達成していると考える。なお、実用化のための最適なパラメータ設定は、実際の計測機器の仕様や測定対象の性質に依存して決定されるべきものであるため、これらの決定のためには実際の計測機器の仕様や測定対象の性質を知る必要がある。
【実施例３】
【００７９】
次に、実施例３では、上記実施例１，２と異なる動き追従Ｘ線ＣＴアルゴリズムについて説明する。
（１）確率モデルについて
（１−ａ）尤度
再構成の対象となる物体のＣＴ画像を１次元ベクトルＸ＝[ｘ１，・・・，ｘＪ]（ここで，Ｊは再構成画像の画素数）で表現する。 ここでは、時間に応じて変形もしくは移動する物体に対する投影を考えるため、時間をＫ個のフレームＳ（ｋ）, ｋ∈｛１，・・・，Ｋ｝に分割し、ｋ番目のフレームでのＣＴ画像をＸ（ｋ）＝[ｘ１（ｋ），・・・，ｘＪ（ｋ）]で表現する。全フレームのＸ（ｋ）（ｋ∈｛１，・・・，Ｋ｝）を、纏めてＸ＝｛Ｘ（１），・・・，Ｘ（ｋ）｝で表現する。
時刻ｔにおける物体に対する投影による像をＹ（ｔ）＝[ｙ１（ｔ），・・・，ｙ１（ｔ）] (ここで、Ｉは検出器数)とし、全時刻の投影を、纏めてＤ＝｛Ｙ（１），・・・，Ｙ（Ｔ）｝で表現する。時刻ｔがｔ∈Ｓ（ｋ）を満たす場合、ＣＴ画像Ｘ（ｋ）の投影を得たと想定するので、光電変換で生じるノイズが独立なポアソン分布に従うことを考えたとき、下記数式１５，数式１６で表現できる。【００８０】
【数１５】
0006044046.tif 000016
【００８１】
【数１６】
0006044046.tif 000017
【００８２】
上記数式（１６）におけるＬ（ｔ）は、ｍｘｎの行列であり、その（ｉ，ｊ）要素は照射されたＸ線のうち検出器ｉに入射する光線が、再構成画像の画素ｊとどれだけ交差するかを表している。また、上記数式（１６）におけるｂは、ブランクスキャン時のフォトン数(照射されたフォトン数)である。ここで、ノイズの分布にはガウス分布を採用することも可能である。【００８３】
（１−ｂ）事前分布
事前分布を下記数式１７で表現する。【００８４】
【数１７】
0006044046.tif 000018
【００８５】
ここで、Ｚ，Ｗは各フレームのＣＴ画像の画素それぞれに対応する変数からなり、それぞれ組織クラスを表す変数Ｚ＝｛Ｚ（１），・・・，Ｚ（Ｋ）｝,
分布形状を制御する変数Ｂ＝｛Ｂ（１），・・・，Ｂ（Ｋ）｝である。さらに、ｋ番目のフレームにおける組織クラス, 分布形状を制御する変数は、それぞれＺ（ｋ）＝[ｚ１（ｋ），・・・，ｚＪ（ｋ）]，Ｂ（ｋ）＝[ｂ１（ｋ），・・・，ｂＪ（ｋ）]と表現される。想定する組織クラスの数をＣとしたとき、各画素の組織クラスは、Ｃ次元ベクトルｚｊ（ｋ）＝[ｚｊ１（ｋ），・・・，ｚｊＣ（ｋ）]で表現される。ｚｊＣ（ｋ）は、そのｋ番目のフレームの画素ｊがとる組織クラスｃに対してのみｚｊＣ（ｋ）＝１をとり、その他の要素はゼロをとる。
したがって、下記数式１８の制約を満足する。また、ｋ番目のフレームの画素ｊの分布形状を制御する変数ｂｊ（ｋ）は、は、下記数式１９の２次元ベクトルとしている。なお、より複雑な分布を表現するためには、ｄの取りうる値の数を２より大きな任意の数に増やすことで対応可能である。【００８６】
【数１８】
0006044046.tif 000019
【００８７】
【数１９】
0006044046.tif 000020
【００８８】
パラメータθは、動きを表すパラメータとして定義され、事前分布は、下記数式２０で表現される。
数式２０におけるＥ（Ｘ，Ｚ，Ｂ，θ）は、下記数式２１で表現される。【００８９】
【数２０】
0006044046.tif 000021
【００９０】
【数２１】
0006044046.tif 000022
【００９１】
以下に、パラメータおよび変数の説明をする。
まず、下記数式２２で表現されるパラメータは、ｚｊｃ（ｋ）＝ｚｓｃ´（ｋ）＝ｂｊｄ（ｋ）＝ｂｓｄ´（ｋ）＝１のときのみ、値Ｇｘ，ｃｄｓｐａｔｉａｌをとる。また、下記数式２３で表現されるパラメータは、ｚｊｃ（ｋ）＝ｚｓｃ´（ｋ）＝１のときのみ、値Ｇｘ，ｃｓｐａｔｉａｌをとる。【００９２】
【数２２】
0006044046.tif 000023
【００９３】
【数２３】
0006044046.tif 000024
【００９４】
したがって、Ｇｘ，ｃｄｓｐａｔｉａｌを変更することによって、隣接画素の組織クラスと分布形状制御変数が同一だった時に、画素値ｘがどれだけ滑らかとなるべきかを制御できる。また、Ｇｘ，ｃｓｐａｔｉａｌを変更することによって、分布形状制御変数によらず隣接画素の組織クラスが同一だった時に画素値xがどれだけ滑らかとなるべきかを制御できる。同様にＧｚ，ｃｓｐａｔｉａｌは隣接画素がどれだけ同じ組織クラスとなりやすいか、Ｇｂ，ｃｄｓｐａｔｉａｌは隣接画素が同じ組織クラスだった時にどれだけ同じ分布形状制御変数をとるかを制御する。Ｇｃｄｓｅｌｆは組織クラス変数、分布形状制御変数のどのような組み合わせが生じやすいかを表す。Ｇｘｔｅｍｐｏｒａｌ, Ｇｚｔｅｍｐｏｒａｌはそれぞれ変数ｘ、ｚがどれだけ動きパラメータθに沿って時間的に遷移しやすいかを制御する変数である。η（ｊ），Ｎ（ｊ）は、画素ｊの近傍の画素インデックスの集合を表す。それぞれ滑らかさをどれだけの範囲に及ぼすかを決定するものであり、計算量の節約のためには小さくとることが好ましい。【００９５】
また、関数Ｇ（ｉ，ｊ，ｍ（ｋ），（θ））は、下記数式２４で定義される。ここで、ｒｉは画素ｉの位置ベクトル、ｍｉ（ｋ）（θ）は画素ｉのフレームｋからｋ＋１への移動量を表す動きベクトルである。つまり、ｋ番目のフレームの画素ｉは（ｋ＋１）番目のフレームにおいて、位置ｒｉから位置ｒｉ＋ｍｉ（ｋ）（θ）に移動すると期待されることを表現する。なお、γはある正の定数である。したがって、関数Ｇ（ｉ，ｊ，ｍ（ｋ），（θ））は、期待される移動量に近い画素ほど関数Ｇ（ｉ，ｊ，ｍ（ｋ），（θ））が大きな値をとるような重み付け関数として働くことになる。【００９６】
【数２４】
0006044046.tif 000025
【００９７】
ここで、移動量を表す動きベクトルｍｉ（ｋ）は、下記数式２５で定義される。これは、ｅｘｐ（−α||ｒｉ−ｒｊ||−β｜ｋ−ｌ｜）が、時空間的な重みと考えれば、動きベクトルｍｉ（ｋ）がパラメータθｉ（ｌ）の補間で得られることを表している。α,βはそれぞれ空間的補間と時間的補間の強さを制御するパラメータである。【００９８】
【数２５】
0006044046.tif 000026
【００９９】
下記数式２６は、Ｇｘ，ｃｄｓｐａｔｉａｌ＝Ｇｘ，ｃｓｐａｔｉａｌ＝Ｇｘｔｅｍｐｏｒａｌ＝０が成り立つ時、下記数式２７のように、Ｐ（Ｘ｜Ｚ，Ｂ，θ）が、Ｚ，Ｂが与えられた下で独立という条件付き独立な分布として表現される。このとき、Ｐ（Ｘ｜Ｚ，Ｂ，θ）は数式２７を改めて下記数式２８のように表現できる。【０１００】
【数２６】
0006044046.tif 000027
【０１０１】
【数２７】
0006044046.tif 000028
【０１０２】
【数２８】
0006044046.tif 000029
【０１０３】
また、上記数式２７におけるｐ（Ｚ，Ｂ｜θ）は下記数式２９で、上記数式２６におけるＥ（Ｘ，Ｚ，Ｂ，θ）は下記数式３０で定義される。このように定義することで設定すべきパラメータを減らすことができる。【０１０４】
【数２９】
0006044046.tif 000030
【０１０５】
【数３０】
0006044046.tif 000031
【０１０６】
Ｐ（ｘｊ（ｋ）｜ｚｊｃ（ｋ）＝１，ｂｊｄ（ｋ）＝１，θ）に関しては、ＣＴ値ｘが非負の値をとることを考慮して、ガンマ分布を利用することも可能である。このとき、Ｐ（ｘｊ（ｋ）｜ｚｊｃ（ｋ）＝１，ｂｊｄ（ｋ）＝１，θ）は、下記数式３１で表現される。数式３１において、ｕｃｄ，ｖｃｄはそれぞれガンマ分布のパラメータである。【０１０７】
【数３１】
0006044046.tif 000032
【０１０８】
次に、上記の確率モデルを用いて、動きのパラメータθを推定しつつ、ＣＴ画像Ｘ, 組織クラス変数Ｚ,  分布形状制御変数Ｂを推定する方法について説明する。
（２）自由エネルギー最小化に基づく推定方法について
表記を簡単にするために観測変数ｙから隠れ変数ｘを推定する一般的な問題を考える。上記の確率モデルでは、観測変数をＤ，隠れ変数をＸ，Ｚ，Ｂと考えれば対応がつく。また、自由エネルギーＦは、下記数式３２で定義される。数式３２において、ｑ（ｘ）は試験分布と呼ばれる分布である。【０１０９】
【数３２】
0006044046.tif 000033
【０１１０】
対数尤度ｌｏｇｐ（ｙ｜θ）は、自由エネルギーＦを用いると下記数式３３で表現できる。ここで、ＫＬ[ｑ（ｘ）｜ｐ（ｘ）]は、ｘが実数の場合は数式３４で定義される関数であり、ｘが離散変数の場合は数式３５で定義される関数であり、確率分布間の近さを測る尺度として用いられる。これは常に非負の値をとり、ｑ（ｘ）＝ｐ（ｘ）が成り立つときのみゼロをとる。【０１１１】
【数３３】
0006044046.tif 000034
【０１１２】
【数３４】
0006044046.tif 000035
【０１１３】
【数３５】
0006044046.tif 000036
【０１１４】
また、パラメータ推定では、しばしば最尤推定が用いられる。最尤推定の目的は下記数式３６の最適化である。上記の式変形からわかるように対数尤度は自由エネルギーを用いて、下記数式３７と表現できる。ここで、数式３７における左辺は、試験分布ｑ（ｘ）に依存しない量であるので、右辺も試験分布のとり方に依存しない。【０１１５】
【数３６】
0006044046.tif 000037
【０１１６】
【数３７】
0006044046.tif 000038
【０１１７】
また、数式３８であるから、数式３９が成り立つ。この式から右辺は、ＫＬ[ｑ（ｘ）｜ｐ（ｘ｜ｙ，θ）]＝０のとき、すなわち、ｑ（ｘ）＝ｐ（ｘ｜ｙ，θ）のとき、最小値−ｌｏｇｐ（ｙ｜θ）が達成されることがわかる。つまり、対数尤度を最大化するパラメータは下記数式４０と表現できる。【０１１８】
【数３８】
0006044046.tif 000039
【０１１９】
【数３９】
0006044046.tif 000040
【０１２０】
【数４０】
0006044046.tif 000041
【０１２１】
したがって、左辺の自由エネルギーを試験分布ｑ（ｘ）に関して最小化したうえで、パラメータθに関して最小化することで最尤推定が行えることがわかる。この事実より、対数尤度の最大化のかわりに自由エネルギーを試験分布ｑ（ｘ）とパラメータθについて最大化するという方法が考えられる。
ここでは、数式４１を求めるために、数式４２に示す２つのステップを繰り返す。
試験分布ｑ（ｘ）の最適値は、事後分布ｐ（ｘ｜ｙ，θ）で与えられるが、一般に事後分布の解析的な計算は困難である。そこで、ｑ（ｘ）の候補を制約することで最適化を可能にする。【０１２２】
【数４１】
0006044046.tif 000042
【０１２３】
【数４２】
0006044046.tif 000043
【０１２４】
ここでは、下記数式４３なる独立性を仮定した。この仮定のもとで、交互最適化によって最適化することを考えると、下記数式４４になる。
そして、変分法によって、数式４５の解を解析的に求めることができる。一方、数式４６に関しては解析的に求まらないので、ｑ（ｘｊ（ｋ））として、下記数式４７のガウス分布を仮定するか、或いは、下記数式４８のガンマ分布を仮定する。【０１２５】
【数４３】
0006044046.tif 000044
【０１２６】
【数４４】
0006044046.tif 000045
【０１２７】
【数４５】
0006044046.tif 000046
【０１２８】
【数４６】
0006044046.tif 000047
【０１２９】
【数４７】
0006044046.tif 000048
【０１３０】
【数４８】
0006044046.tif 000049
【０１３１】
上記のようにすることで、数式４９を解析的に評価することができる。また、数式４９において、ｑ（Ｘ）のパラメータ（ｑ（ｘｊ（ｋ））が、ガウス分布に従うとすれば平均μと分散σ２、ガンマ分布に従うとすればパラメータ（ρｃｄ，ωｃｄ）に関して共役勾配法などの非線形最適化法で最適化することができる。なお、ガンマ分布はｘの非負性を表現できる好ましい性質を有する。【０１３２】
【数４９】
0006044046.tif 000050
【産業上の利用可能性】
【０１３３】
本発明は、医用ならびに産業用のＸ線ＣＴ装置に有用である。医用画像とくに心臓のＸ線ＣＴが重要な応用であるが、産業用の非破壊検査においても微小な物体の測定時には、ステージの微小な動きが問題になったり、流体が存在する場合にはその流体の変形が問題になったりすることから、様々な用途において、本発明のモーションアーティファクトの低減化が役立つことが期待される。
【図１】
0006044046.tif 000051
【図２】
0006044046.tif 000052
【図３】
0006044046.tif 000053
【図４】
0006044046.tif 000054
【図５】
0006044046.tif 000055
【図６】
0006044046.tif 000056
【図７】
0006044046.tif 000057
【図８】
0006044046.tif 000058
【図９】
0006044046.tif 000059
【図１０】
0006044046.tif 000060
【図１１】
0006044046.tif 000061
【図１２】
0006044046.tif 000062
【図１３】
0006044046.tif 000063
【図１４】
0006044046.tif 000064
【図１５】
0006044046.tif 000065