(57)【要約】【課題】高速に収束しかつ精確な結果を得ることの可能な画像再構成法およびこの方法を実行するコンピュータプログラム命令を格納するコンピュータ読み出し可能媒体ならびに画像再構成装置を提供すること。【解決手段】サイノグラムから画像を再構成する方法において、(ａ）  現サイノグラム残余の減アーチファクト再構成画像を画像空間にて形成し、(ｂ）  この現サイノグラム残余と、この現サイノグラム残余の減アーチファクト再構成画像のラドン変換との差分をラドン空間にて最小化して次サイノグラム残余を計算し、（ｃ）  この次サイノグラム残余が収束する場合、ステップ(ａ）で形成した、少なくとも１つの減アーチファクト再構成画像から最終的な再構成結果を形成することによって解決される。【選択図】図１
2007244871.tif 000002
【特許請求の範囲】
【請求項１】
サイノグラムから画像を再構成する方法において、
(ａ）  現サイノグラム残余の減アーチファクト再構成画像を画像空間にて形成し、
(ｂ）  該現サイノグラム残余と、当該現サイノグラム残余の減アーチファクト再構成画像のラドン変換との差分をラドン空間にて最小化して次サイノグラム残余を計算し、
（ｃ）  該次サイノグラム残余が収束する場合、ステップ(ａ）で形成した、少なくとも１つの減アーチファクト再構成画像から最終的な再構成結果を形成することを特徴とする、
サイノグラムから画像を再構成する方法。
【請求項２】
さらに前記の次サイノグラム残余が収束するまで次サイノグラム残余を現サイノグラム残余として使用してステップ（ａ）および（ｂ）を反復して繰り返す、
請求項1に記載の方法。
【請求項３】
ステップ（ｃ）には、ステップ（ａ）の各反復にて形成される減アーチファクト画像を結合することが含まれる、
請求項２に記載の方法。
【請求項４】
ステップ（ａ）および（ｂ）の最初の反復に対して前記現サイノグラム残余はサイノグラム全体である、
請求項２に記載の方法。
【請求項５】
ステップ（ａ）には、画像空間の減アーチファクト再構成画像に事前知識の制約を強制することが含まれる、
請求項1に記載の方法。
【請求項６】
ステップ（ａ）には、
フィルタ逆投影を使用して現サイノグラム残余の画像を再構成し、
フィルタ逆投影を使用して再構成された当該画像を適用的にバイラテラルフィルタリングして事前の制約に基づいて当該画像を正則化することにより、現サイノグラム残余の減アーチファクト再構成画像を形成する、
請求項1に記載の方法。
【請求項７】
ステップ（ｂ）には、
前記の現サイノグラム残余の減アーチファクト再構成画像のラドン変換を計算することによって更新方向を決定し、
当該更新方向に基づいてステップサイズを決定して、前記の現サイノグラム残余と、当該の現サイノグラム残余の減アーチファクト再構成画像のラドン変換との間の差分を最小化し、
前記の現サイノグラム残余から得られるステップサイズと更新方向との積を減算することによって次サイノグラム残余を計算することが含まれる、
請求項1に記載の方法。
【請求項８】
ステップ（ｃ）には、前記の現残余に対して決定されるステップサイズによって、現サイノグラム残余の減アーチファクト再構成画像を重み付けすることが含まれる、
請求項７に記載の方法。
【請求項９】
前記の現サイノグラム残余の減アーチファクト再構成画像は、画像空間における更新方向である、
請求項１に記載の方法。
【請求項１０】
前記のサイノグラムには、コンピュータトモグラフィスキャンから得られた投影測定値が含まれる、
請求項１に記載の方法。
【請求項１１】
サイノグラムから画像を再構成する方法を実行するためのコンピュータプログラム命令を格納するコンピュータ読出可能媒体において、
前記コンピュータプログラム命令を定めているステップは、
(ａ）  現サイノグラム残余の減アーチファクト再構成画像を画像空間にて形成し、
(ｂ）  該現サイノグラム残余と、当該現サイノグラム残余の減アーチファクト再構成画像のラドン変換との差分をラドン空間にて最小化して次サイノグラム残余を計算し、
（ｃ）  該次サイノグラム残余が収束する場合、ステップ(ａ）で形成した、少なくとも１つの減アーチファクト再構成画像から最終的な再構成結果を形成することとを含んでいることを特徴とする、
サイノグラムから画像を再構成する方法を実行するためのコンピュータプログラム命令を格納するコンピュータ読出可能媒体。
【請求項１２】
さらに前記の次サイノグラム残余が収束するまで次サイノグラム残余を現サイノグラム残余として使用してステップ（ａ）および（ｂ）を反復して繰り返すステップを定めるコンピュータプログラム命令を含んでいる、
請求項1１に記載のコンピュータ読み出し可能媒体。
【請求項１３】
ステップ（ｃ）を定めているコンピュータプログラム命令には、ステップ（ａ）の各反復において形成される減アーチファクト画像を結合するステップを定めるコンピュータプログラム命令が含まれる、
請求項１２に記載のコンピュータ読み出し可能媒体。
【請求項１４】
ステップ（ａ）および（ｂ）の最初の反復に対して前記現サイノグラム残余はサイノグラム全体である、
請求項１２に記載のコンピュータ読み出し可能媒体。
【請求項１５】
ステップ（ａ）を定めているコンピュータプログラム命令には、
画像空間の減アーチファクト再構成画像に事前知識の制約を強制するステップを定めるコンピュータプログラ命令が含まれている、
請求項１１に記載のコンピュータ読み出し可能媒体。
【請求項１６】
前記のコンピュータプログラム命令を定めているステップ（ａ）には、
フィルタ逆投影を使用して現サイノグラム残余の画像を再構成し、
フィルタ逆投影を使用して再構成された当該画像を適用的にバイラテラルフィルタリングして事前の制約に基づいて当該画像を正則化することにより、現サイノグラム残余の減アーチファクト再構成画像を形成するステップを定めるコンピュータプログラム命令が含まれている、
請求項１1に記載のコンピュータ読み出し可能媒体。
【請求項１７】
前記のコンピュータプログラム命令を定めているステップ（ｂ）には、
現サイノグラム残余の減アーチファクト再構成画像のラドン変換を計算することによって更新方向を決定し、
当該更新方向に基づいてステップサイズを決定して、前記の現サイノグラム残余と、当該の現サイノグラム残余の減アーチファクト再構成画像のラドン変換との間の差分を最小化し、
前記の現サイノグラム残余から得られるステップサイズと更新方向との積を減算することによって次サイノグラム残余を計算するステップを定めるコンピュータプログラム命令が含まれる、
請求項1１に記載のコンピュータ読み出し可能媒体。
【請求項１８】
ステップ（ｃ）を定めているコンピュータプログラム命令には、前記の現残余に対して決定されるステップサイズによって、現サイノグラム残余の減アーチファクト再構成画像の重み付けを行うステップを定めるコンピュータプログラム命令が含まれる、
請求項１７に記載のコンピュータ読み出し可能媒体。
【請求項１９】
サイノグラムから画像を再構成する装置において、
現サイノグラム残余の減アーチファクト再構成画像を画像空間にて形成する手段と、
該現サイノグラム残余と、当該現サイノグラム残余の減アーチファクト再構成画像のラドン変換との差分をラドン空間にて最小化して次サイノグラム残余計算する手段と、
該次サイノグラム残余が収束する場合、前記の少なくとも１つの減アーチファクト再構成画像から再構成結果を形成する手段とを有することを特徴とする、
サイノグラムから画像を再構成する装置。
【請求項２０】
前記の現サイノグラム残余は、サイノグラム全体である、
請求項１９に記載の装置。
【請求項２１】
前記の再構成結果を形成する手段は、
複数のサイノグラムから形成される減アーチファクト再構成画像に対する手段を含んでいる、
請求項１９に記載の装置。
【請求項２２】
前記の現サイノグラム残余の減アーチファクト再構成画像を画像空間に形成する手段は、
画像空間の当該減アーチファクト再構成画像に事前知識の制約を強制する手段を含んでいる、
請求項１９に記載の装置。
【請求項２３】
前記の現サイノグラム残余の減アーチファクト再構成画像を画像空間に形成する手段は、
フィルタ逆投影を使用して当該現サイノグラム残余の画像を再構成する手段と、
フィルタ逆投影を使用して再構成された当該画像を適用的にバイラテラルフィルタリングして、事前の制約に基づいて当該画像を正則化し、これによって現サイノグラム残余の減アーチファクト再構成画像を形成する手段とを含んでいる、
請求項１９に記載の装置。
【請求項２４】
前記の現サイノグラム残余と、当該現サイノグラム残余の減アーチファクト再構成画像との差分をラドン空間にて最小化する手段は、
前記の現サイノグラム残余の減アーチファクト再構成画像のラドン変換を計算してラドン空間における更新方向を決定する手段と、
当該更新方向に基づいてステップサイズを決定して、前記の現サイノグラム残余と、当該の現サイノグラム残余の減アーチファクト再構成画像のラドン変換との間の差分を最小化する手段と、
前記の現サイノグラム残余から得られるステップサイズと更新方向との積を減算することによって次サイノグラム残余を計算する手段とを含んでいる、
請求項１９に記載の装置。
【請求項２５】
前記の再構成結果を形成する手段は、
前記の現サイノグラム残余の減アーチファクト再構成画像をステップサイズで重み付けする手段と、
当該の乗算結果を前の再構成結果に加える手段とを含んでいる、
請求項２４に記載の装置。
【請求項２６】
前記の現サイノグラム残余の減アーチファクト再構成画像を画像空間に形成する手段は、
当該画像空間における更新方向を決定する手段を含んでいる、
請求項１９に記載の装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、請求項１の上位概念に記載された、サイノグラムから画像を再構成する方法、請求項１１の上位概念に記載されたこの方法を実行するコンピュータプログラム命令を格納するコンピュータ読み出し可能媒体、ならびに請求項１９の上位概念に記載された、サイノグラムから画像を再構成する装置に関する。本出願は、２００６年３月１６日に提出した米国暫定特許明細書第６０／７８３，１１７７号に優先権を主張するものであり、その開示内容は参照によって本願に含まれるものとする。
【背景技術】
【０００２】
本発明は、画像再構成に関し、また殊に医用イメージングにおけるトモグラフィ画像再構成に関する。
【０００３】
トモグラフィ画像再構成とは、ラドン変換によって表される投影測定値の集合から画像を再構成することである。ラドン変換により、測定した投影データからトモグラフィ画像を再構成するための数学的なベースが得られる。画像のラドン変換は、ふつうサイノグラムと称される。トモグラフィ画像再構成では、投影データのサイノグラムに基づいてピクセル値が決定されることによって画像が再構成される。このことはA. KakおよびM. SlaneyによるPrinciples of Computerized Tomographic Images, IEEE Press, 1988に極めて詳細に説明されている。
【０００４】
実際には投影測定値の個数および質は制限されていることが多く、このため、これらの投影測定値から再構成される画像の品質が影響を受ける。例えば、医用イメージングでは、患者への過剰な放射のリスクを低減するため、低線量の診断トモグラフィイメージングが広く使用されている。このため投影測定値の個数および質が大きく制限され得るのである。ふつう非反復形の再構成アルゴリズムにより、上記のような制限された投影測定値が十分に扱われることはなく、重大な縞模様状のアーチファクトが形成されることが多い。このような非反復形アルゴリズムの１つは、広く知られているフィルタ逆投影（ＦＢＰ = Filtered Back Projection）アルゴリズムであり、これは、A. KakおよびM. SlaneyによるPrinciples of Computerized Tomographic Images, IEEE Press, 1988に記載されている。そこでは非反復形アルゴリズムにおけるアーチファクトを低減するためにラドン空間における補間が提案されている。しかしながら投影の零空間において画像形状および境界の事前知識を適用することはふつう極めて困難である。
【０００５】
反復法は、画像空間における事前知識を組み込んで反復処理方式における再構成を正則化することによって非反復法よりもよい結果を達成している。このような反復法においてふつう目的関数は、測定したサイノグラムと再構成画像の前投影との間の類似性を強制するデータフィデリティ項（data fidelity term）と、信号についての事前知識（例えば、鋭いエッジを有する滑らかな曲面）を強制する正則化項とから構成される。しかしながらこの目的関数の最適化はふつう勾配降下法をベースにしている。データフィデリティ項の結果に対する微分は、各反復における逆投影であり、これによってふつうはぼやけた結果しか得られず、また収束も遅い。さらに、勾配降下方向における更新を行うためのステップサイズは、収束の問題のためにふつう小さく設定され、これによって収束に必要な反復数がさらに大きく増大してしまうのである。
【非特許文献１】A. Kak, M. Slane, Principles of Computerized Tomographic Images, IEEE Press, 1988
【発明の開示】
【発明が解決しようとする課題】
【０００６】
本発明の課題は、高速に収束しかつ精確な結果を得ることの可能な画像再構成法およびこの方法を実行するコンピュータプログラム命令を格納するコンピュータ読み出し可能媒体ならびに画像再構成装置を提供することである。
【課題を解決するための手段】
【０００７】
上記の画像再構成法についての課題は、本発明の請求項１により、サイノグラムから画像を再構成する方法において、(ａ）  現サイノグラム残余の減アーチファクト再構成画像を画像空間にて形成し、(ｂ）  この現サイノグラム残余と、この現サイノグラム残余の減アーチファクト再構成画像のラドン変換との差分をラドン空間にて最小化して次サイノグラム残余を計算し、（ｃ）  この次サイノグラム残余が収束する場合、ステップ(ａ）で形成した、少なくとも１つの減アーチファクト再構成画像から最終的な再構成結果を形成することによって解決される。
【０００８】
上記のコンピュータ読み出し可能媒体についての課題は、本発明の請求項１１により、サイノグラムから画像を再構成する方法を実行するためのコンピュータプログラム命令を格納するコンピュータ読出可能媒体において、上記のコンピュータプログラム命令を定めているステップは、(ａ）  現サイノグラム残余の減アーチファクト再構成画像を画像空間にて形成し、(ｂ）  この現サイノグラム残余と、この現サイノグラム残余の減アーチファクト再構成画像のラドン変換との差分をラドン空間にて最小化して次サイノグラム残余を計算し、（ｃ）  この次サイノグラム残余が収束する場合、ステップ(ａ）で形成した、少なくとも１つの減アーチファクト再構成画像から最終的な再構成結果を形成することとを含んでいることを特徴とする、コンピュータ読出可能媒体によって解決される。
【０００９】
上記の画像再構成装置についての課題は、本発明の請求項１９により、サイノグラムから画像を再構成する装置において、現サイノグラム残余の減アーチファクト再構成画像を画像空間にて形成する手段と、この現サイノグラム残余と、この現サイノグラム残余の減アーチファクト再構成画像のラドン変換との差分をラドン空間にて最小化して次サイノグラム残余計算する手段と、この次サイノグラム残余が収束する場合、上記の少なくとも１つの減アーチファクト再構成画像から再構成結果を形成する手段とを有することを特徴とする、サイノグラムから画像を再構成する装置によって解決される。
【発明を実施するための最良の形態】
【００１０】
本発明によって提供されるのは、高速に収束しかつ精確な結果を提供する画像再構成法および画像再構成システムである。本発明の実施形態では、反復形画像再構成が２つのステージ、すなわち画像空間における１ステージと、ラドン空間における１ステージとに分割される。
【００１１】
本発明の１実施形態では、画像空間において適応形フィルタリングを行って、現サイノグラム残余（current sinogram residue）の減アーチファクト再構成画像（artifact-reduced recontructed image）、すなわちアーチファクトを低減した再構成画像を形成する。現サイノグラム残余のこの減アーチファクト再構成画像は、画像空間における更新または降下方向を表す。上記の適応形フィルタリングによって、現サイノグラム残余の再構成画像が、先行の制約に基づいて正則化される。上記の更新方向は、ラドン空間に変換され、ステップサイズが決定される。ここでこのステップサイズは、現サイノグラム残余と、この現サイノグラム残余の減アーチファクト再構成画像のラドン変換との間の差分をラドン空間において最小化するステップサイズである。現サイノグラム残余と、この現サイノグラム残余の減アーチファクト再構成画像のラドン変換との間の最小化されたこの差分が、次サイノグラム残余である。これらのステップは、解が収束するまで次サイノグラム残余毎に反復して繰り返すことができる。最終的な再構成結果（画像）は、サイノグラム残余毎に形成した減アーチファクト再構成画像を結合することによって形成される。
【００１２】
本発明の上記および別の利点は、以下の詳細な説明および添付の図面を参照することによって当業者には明らかになろう。
【実施例】
【００１３】
本発明は、トモグラフィ画像再構成法に関する。本発明は、コンピュータトモグラフィ（ＣＴ）、磁気共鳴イメージング（ＭＲＩ）などを含む任意の医用イメージング法に適用可能である。本発明の複数の実施形態をここに説明して、本発明の信号ダウンサイジング法が視覚的にわかるようにする。これらの実施形態は、コンピュータシステムに格納されたデータを使用して、このコンピュータシステム内で実行できることはいうまでもないことである。したがって、上記の方法のいくつかのステップはコンピュータシステム内の内部表現として発生するものであってもよい。
【００１４】
伝統的な反復形画像再構成法ではラドン変換演算子が離散化されて、各画像ピクセル輝度と、各サイノグラム値との間の重み付けファクタが形成される。重み付けファクタのこの集合は、システムマトリクスＡとして表され、また画像からサイノグラムへの前投影を表す。既知のサイノグラム値が点群をなしてベクトルｂとなることによって、上記の再構成問題は、再構成する未知の画像を表すベクトルｘの線形系を解くこと、すなわちＡｘ＝ｂとして定式化することができる。逆問題が非適切な性質を有することに起因して有効な解を保証するため、滑らかな曲面およびエッジモデリングのような信号特性に基づく正則化がふつう必要である。正則化により、再構成する信号についてのいくつかの事前知識が強制される。したがって目的関数はつぎのように設計することができる。すなわち、
Ｊ ＝ ‖Ａｘ−ｂ‖2＋α‖Ｒｘ‖2                                  （１）
であり、ここで‖Ａｘ−ｂ‖2は、データフィデリティ項、‖Ｒｘ‖2は正則化項、またαはこの正則化項の重み付けである。上記のデータフィデリティ項により、測定されたサイノグラムと、再構成画像の前投影と類似性が強制されて、この再構成画像と真のオブジェクトとが類似することが保証される。しかしながら場合によっては、上記のデータフィデリティ項を満足する複数の解が一般に存在し得る。上記の正則化項によって試みられるのは、画像の事前知識を強制して、この事前知識に適合する解を優先させることである。例えば、一般的な事前知識の強制は、ピクセルがふつう隣接ピクセルと類似の輝度を有することである。このことは、滑らかさについての制約と称することができる。この滑らかさについての制約を強制するため、正則化項を「‖ｘ−近傍における平均輝度‖2」と定めることができる。この項は、滑らかでない解に対して大きなペナルティを与えることになる。
【００１５】
慣用の方法における数式（１）の目的関数の最適化はふつう、微分
∂Ｊ／２∂ｘ ＝ (ＡTＡ＋αＲTＲ)ｘ−ＡTｂ                      （２）
に基づく勾配降下法によって行われ、ここでＡTは逆投影演算子を表し、この演算子によってサイノグラムが画像に再構成される。しかしながら上記の微分から導出される逆投影は、画像の周波数の高い構造を効果的に取り出すことはできず、収束が遅いという欠点を有する。
【００１６】
この最適化のフレームワークにおいてわかるのは、一方では再構成画像を正則化項に基づいて更新して、さまざまな事前知識および制約をこの画像に強制しなければならないことである。これらの事前知識および制約はふつう、ラドン空間よりも画像空間においてより一層効果的に強制される。この一方で複数のサイノグラムを整合させて、上記の目的関数のデータフィデリティ項を強制しなければならないが、これはラドン空間において一層効果的に行い得るのである。慣用の方法は、目的関数のこれらの２つの部分を同時に扱っており、このために最適化が困難なのである。
【００１７】
本発明の１実施形態によれば、トモグラフィ画像再構成は２ステージの方法で実行することができ、この方法では１ステージを画像空間において実行し、また１ステージをラドン空間において実行する。図１には本発明の１実施形態による画像再構成法が示されている。図１に示したようにこの方法は２つのステージ１００および１１０に分割することができる。ステージ１００は、画像空間における適応形フィルタリングステージであり、このステージにより、各反復において再構成画像が改善される。ステージ１２０は、ステージ１００において決定された画像の前投影と、サイノグラム尺度との間の残差を最小化する、ラドン空間における最適化ステージである。図２は、図１の方法の１反復におけるデータの流れを説明するブロック図である。この図を図１に関連して説明する。
【００１８】
ステップ１０２において、ｉ番目の反復に対し、サイノグラム残余sino(i)から画像Ｉ(i)AR２０４を再構成する。初期のサイノグラム残余sino(0)は元々のサイノグラムである。この元々のサイノグラムは、ＣＴスキャナのようなスキャンニング装置によって測定可能な投影測定値である。画像I(i)ARは、サイノグラム残余sino(i)のフィルタ逆投影（ＦＢＰ = Filtered Back Projection）に適応形フィルタリングを使用して得られる減アーチファクト画像である。フィルタ逆投影（ＦＢＰ = Filtered Back Projection）により、サイノグラム残余から画像が再構成されるが、低線量のイメージングでは縞状のアーチファクトが発生し得るため、数式（１）の正則化項Ｒにより、鮮明なオブジェクト境界を保持しながら、これらのアーチファクトを適応的にフィルタリングするための良好なモデリングを得る。減アーチファクト画像I(i)ARは、信号部分空間において、改善され実現可能な更新方向を表す。
【００１９】
ＦＢＰは、サイノグラム残余に整合する画像を供給することによって、このサイノグラム残余と画像空間とを関連付けるために使用される。さらに、ＦＢＰによって逆投影よりも高速に高周波成分が取り出されるため、エッジおよび別の領域による関数的な取り扱いが画像空間において一層容易になる。しかしながらＦＢＰ再構成画像は、前の信号モデリング制約（例えば滑らかさについての制約）を満足せず、付加的な正則化によって直接使用することはできない。したがって本発明の１実施形態では、画像空間適応形フィルタリングを使用して上記の正則化項を強制する。ベイズ画像復元に関連する滑らかさの制約に基づいて、この方法により、（ＦＢＰの結果を直接使用する代わりに）画像
2007244871.tif 000003
が復元される。ここでこの画像により、つぎの目的関数が最適化される。この目的関数において添え字ｐは、ピクセルｐにおいて推定される信号を示す。すなわち、
【数１】
2007244871.tif 000004
であり、ここでλi（ｉ ＝ １，２，３）により、各制約の重み付けが行われ、
2007244871.tif 000005
は、ｐの近傍におけるＦＢＰ再構成画像の平均輝度である。
【００２０】
数式（３）の第１項によって強制されるのは、上記の推定された信号と、ＦＢＰ再構成画像と同じように見えるべきであることである。第２項により、上記の信号の性質がモデリングされ、また滑らかな信号が優先される。ｅpは、０と１との間で変化するピクセルベースのエッジモデリング処理を表す。ｅp ＝ １の場合、滑らかさについての制約は不能にされ、鮮明なエッジが保持される。ペナルティλ3を導入するのは、すべてのピクセルをエッジとして扱うことを回避するためである。
初期のＦＢＰ結果から、アーチファクトを有するオブジェクト境界が得られる。ここではゾーベル（sobel）エッジ検出器をチューニングし、さまざまなサイズおよび閾値を相応に調整してオブジェクト構造のエッジを識別する。別の反復の間に、再びこの処理に進む場合、上記の閾値が変化する。それゆえ、適応形の閾値が使用される。（後の反復において復元可能な）弱いオブジェクト構造と比較してアーチファクトは除去するのが難しいため、最初の反復において閾値を手動でチューニングすることができる。これにより、アーチファクトが可能な限り除去されるが、これにはいくつかの弱い構造も除去されてしまうという犠牲を伴う。またここではゾーベル検出のヒストグラムを推定して、ＣＤＦ（Cumulative Density Function）を探し出して百分位数ベースの閾値を計算することができる。後の反復において、反復毎に上記の閾値を手動でチューニングする代わりにこの百分位数を適応形閾値として使用することができる。
【００２１】
数式（３）の第２項において、平均化に対してあらかじめ定めた重み付けを使用して滑らかな信号を得ると、
2007244871.tif 000006
により、エッジの近くでぼけが発生する。バイラテラルフィルタリング（bilateral filtering）を使用して、
2007244871.tif 000007
を得ることができる。それはこれが効率的であり、かつ非反復であるからである。ふつうのローパスフィルタリングによって、近くのピクセル輝度値が平均化され、またこれに伴って遠くのピクセルに対する重み付けが小さくなるが、バイラテラルフィルタによって第２の次元が加わって、放射の類似性が測定される。この次元によって、中央ピクセル輝度に類似する輝度を有するピクセルに大きな重み付けが割り当てられる。空間変数表現において、バイラテラルフィルタ画像
2007244871.tif 000008
は、２つのフィルタリング成分のガウス実現を使用して、入力ＦＢＰ再構成
2007244871.tif 000009
から得ることでできる。すなわち、
【数２】
2007244871.tif 000010
であり、ここでｃは正規化係数であり、σ1によって空間ローパスフィルタリングが制御される。これに対してσ2により、弁別が制御される。滑らかな領域は、その近傍にまたがって類似のピクセル輝度を共有するため、アーチファクトから得られるピクセル輝度差分は、チューニングしたσ2より格段に小さい。これによって数式（４）の第２の指数部分は１に近くなる。したがってバイラテラルフィルタリングによってアーチファクトが除去されて滑らかになるのである。しかしながらフィルタリング処理がエッジのピクセルに到達した場合、滑らかな領域における近接ピクセルの輝度差分は、σ2よりも格段に大きくなって第２の指数部分は０に近くなる。これによって、これらのピクセルからの寄与分は帳消しにされ、輝度選択的なフィルタリングが適合的に得られるのである。このようにしてエッジがぼけることなくアーチファクトが除去される。
【００２２】
再構成されるオブジェクトの境界をモデリングするため、またエッジ検出閾値thに基づいてピクセルベースの重み付け関数ｅpを定めるため、数式（３）の項を組み合わせて、画像空間における滑らかさについての制約を効率的に得ることができる。ここで上記のエッジ検出閾値は、画像空間においてどのピクセル輝度が表されるかによって適合される。エッジを保持するため、ｅpにより、元々のＦＢＰ再構成結果Ｉfbpの構造境界部において一層大きな重み付けが行われる。さらに、アーチファクトを除去するため、ｅpにより、上記のバイラテラルフィルタ再構成Ibilの非境界部においてに一層大きな重み付けが行われる。これによって上記の減アーチファクト画像ＩARは、
【数３】
2007244871.tif 000011
として形成され、ここでｅpは、閾値thからの輝度差分のチューニング可能な関数となるように選択され、また０と１との間に制限されてつぎのようになる。すなわち、
【数４】
2007244871.tif 000012
であり、ここでｋによってこの関数の峻度が調整され、またはＩfbp,pは、画像Ｉfbpにおけるピクセルpの輝度である。
【００２３】
図１および２に戻ると、ステップ１０４において、減アーチファクト画像Ｉ(i)ARのラドン変換を計算して、ラドン空間における更新方向を決定する。図２に示したように、画像Ｉ(i)AR２０４のラドン変換（または前投影）により、ラドン空間における更新方向であるsino(i)AR２０６が得られる。
【００２４】
ステップ１０６において、反復ｉに対するステップサイズsino(i)ARおよびサイノグラム残余sino(i)に基づいてステップサイズａ(i)を決定する。上記のラドン変換の線形性に基づいて、上記の再構成画像のラドン変換はつぎのように分解される。すなわち、
Ａ(Ｉ＋ａ(i)Ｉ(i)AR) ＝ Ａ(Ｉ)+ａ(i)Ａ(Ｉ(i)AR)                      （７）
である。
【００２５】
この画像空間において、再構成結果lは、各反復の後、つぎのように更新することができる。すなわち、
Ｉ ← Ｉ＋ａ(i)Ｉ(i)AR                                          （８）
であり、次サイノグラム残余推定値sino(i+1)はつぎのように計算可能である。すなわち、
sino(i+1) = sino(i)-a(i)sino(i)AR                                  （９）である。
【００２６】
各反復において、ａ(i)はデータフィデリティ項を最小化するように決定しなければならない。すなわち、
【数５】
2007244871.tif 000013
である。これに相応して、ステップサイズａ(i)は、数式（１０）に表されたフィデリティ項の微分をゼロに等しく設定してａ(i)について解くことによって得ることができ、
【数６】
2007244871.tif 000014
である。
【００２７】
ステップ１０８において次サイノグラム残余sino(i+1)は、更新方向sino(i)ARおよびステップサイズａ(i)に基づいて計算される。図２に示したように更新方向sino(i)AR２０６とステップサイズａ(i)との積はサイノグラム残余sino(i)２０２から減算されて、次サイノグラム残余sino(i+1)２０８が計算される。このことは上の数式（９）に表されている。次サイノグラム残余sino(i+1)は、前のサイノグラム残余の再構成画像Ｉ(i)ARから投影されたサイノグラムsino(i)ARと、目下の前サイノグラム残余sino(i)との残差である。
【００２８】
ステップ１１０においてこの解が収束したか否かが決定される。例えば、上記の次サイノグラム残余sino(i+1)と、閾値とを比較することができる。
次サイノグラム残余sino(i+1)が上記の閾値よりも小さい場合、上記の解が収束したとする。次サイノグラム残余sino(i+1)が上記の閾値よりも小さくない場合、上記のサイノグラム残余がまだ最小化されておらず、上記の解はまだ収束していないとする。この解がまだ収束していない場合、この方法はステップ１１２に進む。この解が収束している場合、この解はステップ１１４に進む。
【００２９】
ステップ１１２においてこの方法はつぎの反復に進み、i+1に対してステップ１０２〜１０８を繰り返す。したがってこの方法はステップ１０２に戻り、反復i+1に対してサイノグラム残余sino(i+1)から画像Ｉ(i+1)ARを再構成する。図２に示したようにＦＢＰおよびアーチファクト低減は、サイノグラム残余sino(i+1)２０８に行われ、その結果、再構成画像Ｉ(i+1)AR２１０が得られる。ステップ１０２〜１０８は、サイノグラム残余と、再構成画像から得られかつ投影されたサイノグラム残余との間の残差が最小化されて解が収束するまで反復して繰り返される。サイノグラム残余は反復毎に減少するため、解が収束することは保証される。
【００３０】
ステップ１１４において、反復毎に形成されるサイノグラム残余の再構成画像を結合することによって最終的な再構成画像が形成される。したがって最終的な再構成結果は、つぎのように表すことができる。すなわち、
Ｉ ＝ ａ(0)Ｉ(0)AR＋Λ＋ａ(i)Ｉ(i)AR＋Λ                        （１２）
である。
したがって最終的な再構成画像は、画像空間における「ステップ」の和であり、ここで各ステップは、降下形の方向ＩARとステップサイズａの積である。
【００３１】
図３（ａ）〜３（ｂ）は、本発明の１実施形態をテストするために使用される投影データの例を示している。図３（ａ）には、１００個の投影および投影当たり３６７個の検出器を有する臨床低線量ＣＴ  Ｐ３０Ａサイノグラム３０２が示されている。再構成画像のサイズは２５８ピクセル×２５８ピクセルである。図３（ｂ）にはサイノグラム３０２の慣用のＦＢＰ投影３０４が示されている。図３（ｂ）に示されているように慣用のＦＢＰ３０４には縞状のアーチファクトが含まれている。
【００３２】
図４（ａ）〜４（ｃ）は、本発明の１実施形態を使用して図３(ａ）のサイノグラム３０２を再構成した結果の例を慣用の反復法と比較して示している。図４（ａ）は、慣用の反復再構成法を１０回反復した後の画像再構成４０２を示している。図４（b）は、慣用の反復再構成法を5０回反復した後の画像再構成４０4を示している。図４（ｃ）には本発明の１実施形態を使用した画像再構成４０６が、１０回の反復の後に示されている。図４（ａ）〜４（ｃ）からわかるように本発明の実施形態により、計算が格段に節約され、アーチファクトが低減され、また空間分解能が犠牲にされることなく、優れた再構成が得られる。
【００３３】
図５は、本発明の１実施形態と、慣用の反復再構成法との性能比較を説明するグラフ５００である。図５に示されているようにプロット５０２は、本発明の１実施形態を示しており、プロット４０５は慣用の反復法を示している。このグラフ５００は、２つの方法の反復毎にラドン空間における残余サイノグラムのエネルギーを計算することによって形成される。グラフ５００によって示されるのは、本発明の実施形態５０２が、慣用の反復法５０４よりも格段に高速に収束し、またこの実施形態によって、より一層良好な解（すなわち、残余サイノグラムに対してエネルギーが低い）が得られることである。
【００３４】
画像再構成のための上記の方法は、周知のコンピュータプロセッサ、メモリユニット、記憶装置、コンピュータソフトウェアおよび他のコンポーネントを使用するコンピュータ上に実現することができる。このようなコンピュータの高レベルのブロック図が図６に示されている。コンピュータ６０２には、このコンピュータ６０２の全体処理を制御するプロセッサ６０４が含まれており、ここでこの処理は、このような処理を定めるコンピュータプログラム命令を実行することによって行われる。このコンピュータプログラム命令は、記憶装置６１２（例えば磁気ディスク）に記憶することができ、またこのコンピュータプログラム命令の実行が所望される場合にメモリ６１０にロードされる。したがって、画像空間およびラドン空間における上記の方法のステップを実行するアプリケーションは、メモリ６１０および／または記憶装置６１２に格納されるコンピュータプログラム命令によって定めることができ、またこれらのコンピュータプログラム命令を実行するプロセッサ６０４により制御することができる。コンピュータ６０２はまた、ネットワークを介して他の装置と通信するため、１つまたは複数のネットワークインタフェース６０６も有している。コンピュータ６０２にはまた他の入／出力装置６０８が含まれており、これらの入／出力装置により、ユーザとコンピュータ６０２との対話が可能になる（例えばディスプレイ、キーボード、マウス、スピーカ、ボタンなど）。当業者であれば、実際のコンピュータのインプリメンテーションが他のコンポーネントも含み得ること、また図６が、説明のためにこのようなコンピュータのコンポーネントのうちのいくつかを高レベルで表現したものであることを認識できよう。
【００３５】
上記の詳細な説明は、あらゆる点から見て説明および例示であり、制限的なものでないことは無論のことである。またここに開示される本発明の範囲は、上記の詳細な説明によって定められるものではなく、特許法によって許容される最も広い範囲で解釈した請求項によって定められるものである。上で示して説明した実施形態は、本発明の基本原理を説明するためだけのものであり、当業者であれば、本発明の範囲および精神を逸脱することなくさまざまな変更を行い得ることはいうまでもないことである。また当業者であれば本発明の範囲および精神から逸脱することなく他の特徴のさまざまな組み合わせを実現することができよう。
【図面の簡単な説明】
【００３６】
【図１】本発明の１実施形態による画像再構成法を示す図である。
【図２】図１の方法の１反復におけるデータの流れを示すブロック図である。
【図３】本発明の１実施形態をテストするのに使用される投影データの例を示す図である。
【図４】慣用の反復法と比較した本発明の１実施形態の結果の例を示す図である。
【図５】本発明の１実施形態と慣用の反復法との性能比較を示すグラフである。
【図６】本発明を実現可能なコンピュータの高レベルのブロックを示す図である。
【符号の説明】
【００３７】
６０２  コンピュータ、  ６０４  プロセッサ、  ６０６  ネットワークインタフェース、  ６０８  Ｉ／Ｏ、  ６１０  メモリ、  ６１２  記憶装置
【図１】
2007244871.tif 000015
【図２】
2007244871.tif 000016
【図５】
2007244871.tif 000017
【図６】
2007244871.tif 000018
【図３】
2007244871.tif 000019
【図４】
2007244871.tif 000020