(57)【要約】  本発明は、投影測定値miからの対象１の断層画像内の減衰係数μjの反復再構成のための装置及び方法に関する。n番目の反復の間のμjnに対する更新式において、逆投影誤差mi-mi-(μjn)は、ボクセル依存係数式Iにより重み付けされる。このようなボクセルに依存する更新は、特に、ＡＲＴ又はＭＬのアルゴリズムに含まれることができる。
2008532683.tif 000002
【特許請求の範囲】
【請求項１】
投影測定値mi（i=1, ... P）からの断層画像の減衰係数μj（j=1, ... N）の反復再構成のためのデータ処理ユニットにおいて、
所定の測定関数R及び所定の重み行列A=[Aij]を用いて
mi=R(pi)  （１）
及び
【数１】
2008532683.tif 000011
であり、
前記データ処理ユニットが、以下の式
【数２】
2008532683.tif 000012
によってn番目の反復ステップにおいて前記減衰係数μjnの更新値μjn+1（n=1, ...）を決定するように構成される、データ処理ユニット。
【請求項２】
前記測定関数Rが恒等式であることを特徴とする、請求項１に記載のデータ処理ユニット。
【請求項３】
前記インデックセットInが１つの要素のみを含むことを特徴とする、請求項１に記載のデータ処理ユニット。
【請求項４】
前記重み関数Wiが、
【数３】
2008532683.tif 000013
により与えられることを特徴とする、請求項１に記載のデータ処理ユニット。
【請求項５】
前記測定関数Rが、
R(pi)=exp(-pi)  （５）
として定義されることを特徴とする、請求項１に記載のデータ処理ユニット。
【請求項６】
i番目の投影測定値miが、減衰のそれぞれ前及び後の光子の数di及びYiから、
mi=Yi/di  （６）
によって決定されることを特徴とする、請求項１に記載のデータ処理ユニット。
【請求項７】
前記重み関数Wiが、
Wi(A, pin)=exp(pin)  （７）
により与えられることを特徴とする、請求項１に記載のデータ処理ユニット。
【請求項８】
対象の断層画像の生成のための検査装置において、前記対象の投影測定値の生成のためのイメージングデバイスと、請求項１に記載のデータ処理ユニットとを有する検査装置。
【請求項９】
式（１）及び（２）に従う投影測定値mi（i=1, ... P）からの断層画像の減衰係数μj（j=1, ... N）の反復再構成のための方法において、前記減衰係数μjnの更新値μjn+1（n=1, ...）が、式（３）によって決定される、方法。
【請求項１０】
断層画像の減衰係数μj（j=1, ... N）の反復再構成のためのコンピュータプログラムが記憶される記録担体において、前記プログラムが、請求項９に記載の方法を実行するように構成される、記録担体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、投影測定値からの断層画像（tomographic images）の反復再構成（iterative reconstruction）のためのデータ処理ユニット、検査装置、記録担体及び方法に関する。
【背景技術】
【０００２】
対象の断面又は断層画像は、異なる方向から前記対象を通る十分に多数の投影が利用可能である場合に、再構成されることができ、ここで各投影は、投影光線に沿った減衰係数の線積分を表す。たくさんの再構成アルゴリズムが開発され、文献に記載されており、例えば、フィルタ逆投影（Filtered Backprojection）ＦＢＰ型の再構成アルゴリズム（例えばL. A. Feldkamp, L. Davis, and J. Kress: "Practical Cone-beam Algorithm", Journal of Optical Society of America, 1:612-619, 1984）、代数的再構成技術ＡＲＴ（R. Gordon, R. Bender, and G. T. Herman: "Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography", J. Theor. Biol., 29:471-481, 1970）、最尤（ＭＬ）法（K. Lange and J.A. Fessler: "Globally convergent algorithms for maximum a posteriori transmission tomography", IEEE Transactions on Image Processing, 4(10):1430-1450, 1995）、又はこれらの"順序部分集合（Ordered Subsets）"変形例である。これらのアルゴリズムの各々は、特有の長所及び問題点を持つ。したがって、順序部分集合最尤法（ＯＳＭＬ）を用いて再構成された画像は、ＦＢＰ又はＡＲＴを用いる再構成により達成されるより良い信号対雑音比を持つが、ＯＳＭＬを用いる再構成は、ＡＲＴより長い時間及びＦＢＰより大幅に長い時間を必要とする。
【発明の開示】
【発明が解決しようとする課題】
【０００３】
この状況に基づき、本発明の目的は、断層画像の改良された反復再構成に対する手段、特に、高速アクセラレータボード（fast accelerator board）における実施の機会を提供する手段を提供することである。
【課題を解決するための手段】
【０００４】
この目的は、請求項１に記載のデータ処理ユニット、請求項８に記載の検査装置、請求項９に記載の方法及び請求項１０に記載の記録担体により達成される。好適な実施例は、従属請求項に開示される。
【０００５】
第１の態様によると、本発明は、データ処理ユニット、例えば、適切なソフトウェアと一緒に中央処理ユニット、メモリ及びＩ／Ｏインターフェース等を有する従来の（汎用）コンピュータシステムに関する。前記データ処理ユニットは、対象（例えば患者の体）を通って異なる方向から生成された投影測定値mi（i=1, ... P）から前記対象の断層画像の減衰係数μj（j=1, ... N）を反復的に再構成することができるべきである。前記投影測定値は、例えば、ＣＴスキャナにおいてＸ線放射線により生成されることができる。各投影測定値miは、１つの画素及び特定の投影視野に対応する。測定値miは、以下の式、すなわち、
mi=R(pi)  （１）
及び
【数４】
2008532683.tif 000003
によって減衰係数μjに関連し、ここでRは所定の"測定関数"であり、A=[Aij]は所定の重み行列である（[...]が行列を示し、Aijがこの行列の行i及び列jの要素を示し、インデックスi、jの範囲が常に上記のように規定され、すなわち1≦i≦P及び1≦j≦Nであることに注意する）。これら及び以下の式は、図面を参照して下でより詳細に説明される。
【０００６】
前記データ処理ユニットは、n番目の反復ステップの間に以下の式、すなわち、
【数５】
2008532683.tif 000004
によって減衰係数μjnの更新値μjn+1（n=1, ...）を決定するように構成され、ここで変数は以下の意味を持つ、すなわち、
λは緩和パラメータであり、
In⊆{1, ... P}は、例えば"順序部分集合"方法の原理により選択されることができるインデックスセットであり、
Vは、引き数（argument）の単調関数であり、特に恒等式であることができ、
pin=ΣνAiνμνnは、i番目の投影に対する現在の推定値であり、
Wiは、典型的には、基礎をなす再構成モデル（例えばＡＲＴ又はＯＳＭＬ）により決定される所定の重み関数であり、
δin=mi-min=mi-R(pin)は、順投影（forward projection）の現在の誤差である。
【０００７】
上に記載されたデータ処理ユニットは、多くの既知の反復再構成アルゴリズムの基本的なステップ、特に、
１．減衰係数μjnの現在の推定又はモデルに基づく順投影pinの計算と、
２．前記計算された投影と前記測定された投影との間の残留誤差δinの計算と、
３．減衰係数μjnの全域での計算された誤差δinの逆投影又は分配と、
を実施する。
【０００８】
式（３）内の誤差δinの前述の分配は、重み関数Wi、緩和パラメータλ、及び（順序）部分集合Inに対する和のような既知の構成要素を有する。式（３）の重要なフィーチャは、しかしながら、ボクセルに依存する因子V(Aijμjn/Pin)をも有することである。この因子は、各反復ステップにおいて各ボクセルにより経験される更新が、前記ボクセルの対応する順投影に対する寄与に依存することを保証する。シミュレーションが示すように、このようなボクセルに依存する更新は、従来のアルゴリズムの信号体雑音比及び速度を向上する。
【０００９】
測定関数Rは、特に、恒等式であってもよい。この場合、投影測定値miは、典型的にはＡＲＴに基づくアルゴリズムにおいて使用される値piに対応する。
【００１０】
式（３）の一般的定式化は、インデックスセットInが全ての投影を有し、すなわちIn={1, ... P}である場合を含む。真部分集合（genuine subsets）In⊂{1, ... P}は、特に、収束の速度を加速する順序部分集合アプローチにおいて選択される。更に、インデックスセットInは、オプションとして、１つの要素in∈{1, ... P}のみを含んでもよい。n番目の反復ステップに対する１つの（任意の）投影inのこのような選択は、ＡＲＴに対して典型的である。
【００１１】
ＡＲＴに対する式（３）の更なる対応は、重み関数Wiが特に、
Wi(A, Pin)=1/(ΣνAiν2)  （４）
として選択される場合に達成される。
【００１２】
測定関数Rの代替的な定義は、
R(pi)=exp(-pi)  （５）
により与えられる。
【００１３】
Rの前述の定義は、最尤法アプローチに特に適しており、特定の投影幾何構成及び時間において検出器画素に到達する光子の数Yiが測定される。更に、放射線源を出る（すなわち前記対象に入る前の）光子の対応する数diは、典型的には、既知であると仮定される。この場合、投影測定値miは、減衰対象の後ろの光子数Yi及び前記減衰対象の前の光子数diの商により、
mi=Yi/di  （６）
によって定義されることができる。
【００１４】
Yiが、Yi=di・exp(-ΣνAiνμν)として計算されることができ、これが、式（５）が成り立つ場合にYi=di・R(pi)に対応することに注意する。
【００１５】
ＭＬアプローチに対する式（３）の更に他の対応は、重み関数Wiが、
Wi(A, Pin)=exp(pin)  （７）
になるように選択される場合に達成される。
【００１６】
"真の（genuine）"ＭＬ式に対する比較は、式（３）が、より少ない加算を必要とし、したがって計算的により安価であることを示す。
【００１７】
本発明は、対象の断層画像の生成のための検査装置にも関し、前記装置は、上に記載された種類のデータ処理ユニットに結合された前記対象の投影測定値の生成に対するイメージングデバイスを有する。前記イメージングデバイスは、特に、（回転式）Ｘ線デバイス、例えばＣＴスキャナ又はＣアームシステムでありうる。
【００１８】
本発明は、更に、式（１）及び（２）に従う投影測定値mi（i=1, ... P）から断層画像の減衰係数μj（j=1, ... N）の反復再構成のための方法に関し、減衰係数μjnの更新値μjn+1（n=1, ...）は、式（３）によって決定される。前記方法の更なる発展形は、前記データ処理ユニットの好適な実施例に関連して記載された前記ステップの１つ以上を有することができる。したがって、前記方法の詳細、利点及び変形例に関するこれ以上の情報に対して、前記データ処理ユニットの記載が参照される。
【００１９】
最後に、本発明は、記録担体、例えばフロッピー（登録商標）ディスク、ハードディスク又はコンパクトディスク（ＣＤ）を含み、投影測定からの断層画像の減衰係数の反復再構成のためのコンピュータプログラムが前記記録担体に記憶され、前記プログラムは、前述の種類の方法を実行するように構成される。
【００２０】
本発明のこれら及び他の態様は、以下に記載される実施例を参照して説明され、明らかになるだろう。
【００２１】
以下、本発明は、添付の図面を利用して例として記載される。
【発明を実施するための最良の形態】
【００２２】
図１において、断層再構成問題が概略的に図示され、以下に使用される変数が導入される。断層画像再構成の目的は、対象１、例えば患者の体を通るスライス内の減衰係数μの（２又は３次元の）空間分布を決定することである。前記問題の通常の離散的定式化において、前記画像は、インデックスj∈{1, ... N}を持つＮ個のボクセルに分割され、各ボクセルjは、前記減衰係数の対応する値μjを持つ。
【００２３】
放射線源２、例えばＸ線管から、光線iは、対象１を通るように方向付けられる。対象１に吸収されないこれらの光線の光子は、対象１の後ろの検出器３に到達する。検出器３は、典型的には、多数の検出器素子又は画素を有し、前記検出器素子の各々は、投影測定値ｍiを提供する。最尤法（ＭＬ）アプローチにおいて、照射時間の間に画素に到達する光子の数Yiが、典型的には測定され、放射線源２を出る光子の対応する数diは、既知であると仮定される。代数的再構成技術（ＡＲＴ）のようなアプローチにおいて、対数投影値pi=-log(Yi/di)は、典型的には測定値と見なされる。両方のアプローチが、測定変数miにより共通の枠組みで記載されることができ、ここでＭＬアプローチに対してmi=exp(-pi)であり、ＡＲＴアプローチに対してmi=piである。一般に、所定の関数Rを用いてmi=R(pi)である。
【００２４】
１つの投影幾何構成及び１つの投影光線iのみが図１に示されているが、放射線源２及び／又は検出器３が、異なる投影方向を生じるために移動（例えば一般的には回転）され、複数の画素が、前記画素の配置の各々に対して測定されることは言うまでもない。全ての幾何構成及び画素に対する測定されたデータmiは、この場合、インデックスi∈{1, ... P}により線形にインデックスを付けられる。
【００２５】
探される減衰係数μjは、対象１を通る線積分（又は和）として測定値miのみで表され、測定値mi=R(pi)に対する（μjの）各ボクセルjの寄与は、pi=ΣνAiνμνによって重み係数Aijにより与えられる。
【００２６】
行列A=[Aij]は、前記幾何構成及び前記測定の特定のモデルによって事前に計算されることができる（ＡＲＴ又はＭＬに関する文献を参照）。
【００２７】
以下、標準的なＡＲＴ及びＭＬアルゴリズムの短い概説が最初に与えられる。mi=piに対して、測定された投影ベクトルp=[pi]は、行列積
p=Aμ
により表されることができる。
【００２８】
線形連立方程式を解く標準的なＡＲＴ（上記R. Gordon et al.; Y. Censor: "Row-action methods for huge and sparse systems and their applications", SIAM Review, 4(23):444-466, 1981参照）における１つの反復ステップ（n→n+1）は、j番目の吸収係数μj及びi番目の投影piに対して、
【数６】
2008532683.tif 000005
により与えられ、ここでλは、収束速度を制御するのに使用されるいわゆる緩和パラメータであり、δin=(pi-ΣνAiνμνn)=(pi-pin)は、現在の誤差（測定と順投影との間の差）である。
【００２９】
最尤法（ＭＬ）アプローチは、一方、測定される光子のポアソン統計を考慮に入れる。様々な投影iが独立なので、観測された光子カウントYiの対数尤度は、
【数７】
2008532683.tif 000006
と書き表されることができ（上記K. Lange et al.参照）、ここで、diは、i番目の投影に沿った前記線源を出る光子の予期される数であり、Aijは基準体積要素であり（但しΣνAiν=1）、c1は重要でない定数である。μjn>0に対して、前記尤度は、
【数８】
2008532683.tif 000007
の反復ステップ（n→n+1）により増大される。
【００３０】
凸型（convex）ＭＬアルゴリズムの１回の更新は、全ての投影iに対する和を必要とし、これは多くの計算時間がかかり、再構成された画像の遅いが証明されている収束の結果になる。全ての投影を同時に使用するこの要件は、全ての投影の部分集合Inのみを同時に使用することにより軽減されることができる。順序部分集合の選択を用いると、このＭＬ法は、順序部分集合最尤（ＭＬ）法になり、これは、
【数９】
2008532683.tif 000008
と記載されることができ、ここで大域的緩和パラメータλは、収束の速度を制御するために含まれる。
【００３１】
上に記載されたＡＲＴ、ＭＬ又はＯＳＭＬのような通常のアルゴリズムは、n番目の反復（n→n+1）に対する以下の一般的構造、すなわち、
１．前記減衰係数の現在の推定値μnにより与えられる前記画像を通る順投影minを計算（例えば線積分pin又は光子の相対的な数Yin/diを計算）すること、
２．順投影minと測定値miとの間の差δinを計算し、前記モデル（例えばＭＬ又はＡＲＴ）により与えられるように前記差を重み付けすること、
３．前記モデルにより与えられるように前記画像にこれらの重み付けされた差を投影し戻すこと、
を持つ。
【００３２】
本発明によると、専用ハードウェア上で容易に実施されることができる、ステップ３における単純なボクセルに依存する更新が提案される。このような更新の間に、逆投影される差δinは、前記順投影に対する各ボクセルの寄与が後に続く更新に対する重みを設定するように不均一に分配される。このアイデアを使用すると、例えばＡＲＴは、
【数１０】
2008532683.tif 000009
により規定される更新ステップを用いて強化ＡＲＴ（ＥＡＲＴ）に変換されることができ、ここで、（Aijμjn／ΣνAiνμνn）は、前記逆投影に対する各ボクセルの相対的な寄与を設定し、通常のＡＲＴと異ならせる追加の項である。
【００３３】
ＯＳＭＬに対しても、ＯＳＥＭＬと称される強化バージョンが、更新式、
【数１１】
2008532683.tif 000010
を必要とすることにより作成されることができ、ここでexp(-pi)=Yi/diである。
【００３４】
図２は、１回転に対して１１６０の視野を持つ透過断層撮影スキャンをシミュレートするのに使用された１７のビーズ（beads）と、５００ｍｍの直径を持つ（８００のＳＮＲを持つ投影を使用するＦＢＰで再構成された）水ファントムの画像を示す。１つのビーズは、前記ファントムの中心に配置され、８つのビーズが、前記中心から１００ｍｍの距離に配置され、８つのビーズが、前記ファントムの中心から２００ｍｍの距離に配置される。前記ビーズは、スチールのμ値を持つ。前記中心のビーズは、"ビーズ１"と称され、内側の円にある８つのビーズは、"ビーズ２"ないし"ビーズ９"と称され、外側の円の８つのビーズは、"ビーズ１０"ないし"ビーズ１７"と称される。各ビーズの周りのノイズは、ボックスＢ１とボックスＢ２との間のボクセルを使用して計算される。
【００３５】
この空間解像度ファントムを用いて、ＯＳＥＭＬがＯＳＭＬと同じ解像度及びノイズを持ち（図４）、ＥＡＲＴがＡＲＴと同じ解像度において低いノイズを持つ（図３）ことが示されることができる。図３及び４において、測定されるＭＴＦは、（ａ）前記中心のビーズに対して、（ｂ）中心から１００ｍｍ離れて配置された１つのビーズに対して、及び（ｃ）最も外側のビーズの１つに対して示される。また、図（ｄ）は、各ビーズの近隣におけるノイズ（ＳＴＤ）をハウンズフィールド（Hounsfield）単位で表し、これは、前記ファントムの前記中心から同じ距離における全てのビーズに対してほぼ等しい。特に、図３は、ＥＡＲＴ、ＡＲＴ及びＯＳＭＬ再構成に対する結果を表し、図４は、ＯＳＥＭＬ及びＯＳＭＬ再構成に対する結果を表す。両方の場合に、５０のＳＮＲを持つ水解像度ファントムを通る投影が使用された。同様の向上は、８００のＳＮＲを持つ投影を使用して得られることができた。
【００３６】
要約すると、ボクセルに依存する更新を持つ方法が提案され、前記更新は、１つの順投影に沿って不均一に分配される。結果として生じるアルゴリズムは、ＯＳＭＬ再構成の性能を持つ。前記順投影の計算中に、各ボクセルは、この計算に対して異なって寄与する。各ボクセルに対する相対的な重みは、前記順投影により除算された逆投影に対するボクセル寄与により与えられ、前記逆投影中に各ボクセルに対して前記更新を重み付けするのに使用される。前記方法は、ＣＴ及びＸ線再構成において、一般的には、順投影を計算しなければならないあらゆる反復際構成アルゴリズムにおいて使用されることができる。
【００３７】
最後に、本願明細書において、用語"有する"が他の要素又はステップを除外せず、"１つの（"a" or "an"）"が複数を除外せず、単一のプロセッサ又は他のユニットが複数の手段の機能を満たしてもよいことが指摘される。本発明は、一つ一つの新規な特徴的フィーチャ並びに特徴的フィーチャの一つ一つの組み合わせに存在する。更に、請求項内の参照符号は、前記請求項の範囲を制限すると解釈されるべきでない。
【図面の簡単な説明】
【００３８】
【図１】反復再構成問題の一般的な構成を図示する。
【図２】一部の結果が図３及び４に示されるシミュレーション実験に対して使用された１７のスチールビーズを有する水ファントムの断層画像を表す。
【図３ａ】本発明により拡張されたＡＲＴアルゴリズム（ＥＡＲＴ）が標準的なＡＲＴ及びＯＳＭＬアルゴリズムと比較される場合の図２の画像に対する模範的なスチールビーズに対する解像度測定を描く。
【図３ｂ】本発明により拡張されたＡＲＴアルゴリズム（ＥＡＲＴ）が標準的なＡＲＴ及びＯＳＭＬアルゴリズムと比較される場合の図２の画像に対する模範的なスチールビーズに対する解像度測定を描く。
【図３ｃ】本発明により拡張されたＡＲＴアルゴリズム（ＥＡＲＴ）が標準的なＡＲＴ及びＯＳＭＬアルゴリズムと比較される場合の図２の画像に対する模範的なスチールビーズに対する解像度測定を描く。
【図３ｄ】本発明により拡張されたＡＲＴアルゴリズム（ＥＡＲＴ）が標準的なＡＲＴ及びＯＳＭＬアルゴリズムと比較される場合の図２の画像に対する雑音測定を描く。
【図４ａ】本発明により拡張されたＯＳＭＬアルゴリズム（ＯＳＥＭＬ）が標準的なＯＳＭＬアルゴリズムと比較される場合の図３と同じ図を描く。
【図４ｂ】本発明により拡張されたＯＳＭＬアルゴリズム（ＯＳＥＭＬ）が標準的なＯＳＭＬアルゴリズムと比較される場合の図３と同じ図を描く。
【図４ｃ】本発明により拡張されたＯＳＭＬアルゴリズム（ＯＳＥＭＬ）が標準的なＯＳＭＬアルゴリズムと比較される場合の図３と同じ図を描く。
【図４ｄ】本発明により拡張されたＯＳＭＬアルゴリズム（ＯＳＥＭＬ）が標準的なＯＳＭＬアルゴリズムと比較される場合の図３と同じ図を描く。
【図１】
2008532683.tif 000014
【図２】
2008532683.tif 000015
【図３ａ】
2008532683.tif 000016
【図３ｂ】
2008532683.tif 000017
【図３ｃ】
2008532683.tif 000018
【図３ｄ】
2008532683.tif 000019
【図４ａ】
2008532683.tif 000020
【図４ｂ】
2008532683.tif 000021
【図４ｃ】
2008532683.tif 000022
【図４ｄ】
2008532683.tif 000023