JP6101209B2 - コンピュータ断層撮影(ct)のための逆投影の最適化実装 - Google Patents

コンピュータ断層撮影(ct)のための逆投影の最適化実装 Download PDF

Info

Publication number
JP6101209B2
JP6101209B2 JP2013550621A JP2013550621A JP6101209B2 JP 6101209 B2 JP6101209 B2 JP 6101209B2 JP 2013550621 A JP2013550621 A JP 2013550621A JP 2013550621 A JP2013550621 A JP 2013550621A JP 6101209 B2 JP6101209 B2 JP 6101209B2
Authority
JP
Japan
Prior art keywords
weight table
downsampled
weight
computed tomography
perspective geometry
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.)
Active
Application number
JP2013550621A
Other languages
English (en)
Other versions
JP2014506491A (ja
Inventor
アリ, ピー. バーマン,
アリ, ピー. バーマン,
ダニエル ナール,
ダニエル ナール,
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Mercury Systems Inc
Original Assignee
Mercury Systems Inc
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Mercury Systems Inc filed Critical Mercury Systems Inc
Publication of JP2014506491A publication Critical patent/JP2014506491A/ja
Application granted granted Critical
Publication of JP6101209B2 publication Critical patent/JP6101209B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/006Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/421Filtered back projection [FBP]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/428Real-time

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Description

本願は、その内容を引用して本明細書に援用する2012年1月21日付けの米国特許仮出願第61/435,107号および2012年1月20日付けの米国特許仮出願第13/354587号の利益を主張する。
断層像(原語: tomogram)は、三次元物体のスライスすなわち断面の二次元像である。断層写真撮影装置(原語: tomograph)は、断層像を生成するための装置である。断層撮影法(原語: tomography)は、断層像を作成するための工程全体を指す。コンピュータ断層撮影(CT)は、断層像を作成するための工程全体であり、その計算はコンピュータにより実行される。
CTは2つの主たる段階すなわち工程からなる。第1段階は三次元物体の走査であり、第2段階は走査した物体の断層像の計算である。多数の断層像を組み合わせて、走査した物体の三次元密度マップまたは体積表現が作成される。
典型的なCT走査システムは、X線管のような放射線源と、放射線検出器と、コンピュータシステムとから構成されている。放射線源および検出器は、撮像する物体の反対側に配置される。次に、放射線ビームが放射線源から検出器に向かって投射され、物体により吸収されない陽子は、検出器に向かって送出されかつ検出器に衝突する。その結果、X線源の現在位置から三次元物体の二次元投影像を表す像が検出器上に得られる。放射線源および検出器は、典型的には物体の周りを180度または360度回転させつつ、その間に、この撮像処理を複数の中間位置で繰り返すことにより、一定範囲の角度方向にわたる物体の一連の二次元像が得られる。CT走査システムは、物体を囲む一組の固定X線源および検出器を用いて物体を走査することもでき、その場合、X線源も検出器も移動もされない。CT走査は、物体が静止した状態または走査装置に対して直交して運動する状態で行うことができる。後者の場合はヘリカル走査となる。
検出器からの一連のこれら投影像がコンピュータシステムに与えられる。このコンピュータシステムは、これら二次元投影図を用いて物体を様々に再構築できる。この概念は、断層撮影体積再構築と呼ばれている。フェルドカンプ逆投影法、代数的再構築手法(ART)、最尤推定-期待値最大化法(MLEM)を含むがそれらに限定されない様々な数学的アルゴリズムを断層撮影体積再構築と組み合わせて使用できる。ほとんどのアルゴリズムは、多数の投影測定を実行し、物体のすべての点が多くの角度から放射路に含まれることを前提にしている。フェルドカンプ逆投影法は1つの再構築技法であって、投影データがまずフィルタによって畳み込みされ、それぞれのビューは、捕捉の瞬間におけるX線源の角度に対応した角度で、撮像された体積を表す正方形グリッドに連続的に重畳される。重畳または集積の過程において、グリッドの各要素の検出器への投影位置を得るためには透視幾何学的形状(原語:perspective geometry)が既知であることが不可欠であり、乗法重み因子(原語:weight factor)も既知でありかつフィルタされた検出器データからの値に適用されることも不可欠である。
透視幾何学的形状および乗法重み因子は、計算コストが高い超越関数を評価する必要がある。重み因子および透視幾何学的形状は解析関数によって求めることができず、X線源および検出装置の位置の事前情報を用いて求めなければならないので、CT走査システムによってはこれらの計算はさらに複雑となる。いずれの場合も、透視幾何学的形状および重み因子は事前に計算し、メモリに格納しておくことができる。スキャナによっては、検出器およびX線源のアレイは、1つまたは2つの軸に関して完全に対称に構築されているものもある。これらの幾何学的対称性を利用することにより、幾何学的形状および重み因子テーブルのサイズを小さくできる。
物体が移動している場合のCT走査に関しては、スキャナを通る運動が投影の生成と適切に連係されていれば、重み計算もなんらかの周期的間隔で反復される。従って、これにより重みテーブルおよび透視幾何学的形状が、静止物体のCT走査の場合と同様に事前に計算可能となる。
テーブルサイズを縮小するためスキャナの幾何学的形状の対称性を活用したとしても、これらテーブルは非常に大きくなることがあり、メモリストレージおよびメモリアクセス帯域幅要件が再構築システムの性能に悪影響を及ぼす可能性がある。
従って、本発明の目的は、CTにおける逆投影再構築のための事前計算テーブルの格納およびアクセスに関する改良方法を提供することである。この方法により、逆投影法を用いた再構築を実行するCTコンピュータシステムが高速化しかつ/または安価となる。
本発明の実施形態は、所与のコンピュータ断層撮影装置に関連付けられた透視幾何学的形状データおよび重みテーブルをダウンサンプリングしかつ復元するための方法を提供する。前記方法は、前記コンピュータ断層撮影装置に関連付けられた透視幾何学的形状データを生成する段階を含む。前記方法は、前記コンピュータ断層撮影装置に関連付けられた重みテーブルを生成する段階をさらに含む。格納するデータの量を減少させるために、前記重みテーブルおよび前記透視幾何学的形状データがダウンサンプリングされる。前記コンピュータ断層撮影装置に関連付けられた前記ダウンサンプリングされた透視幾何学的形状データおよび前記ダウンサンプリングされた重みテーブルが、記憶装置に格納される。前記方法は、前記コンピュータ断層撮影装置を用いて走査された物体の断層撮影体積の再構築工程の間に、前記透視幾何学的形状データおよび前記重みテーブルを、それぞれ前記ダウンサンプリングされた透視幾何学的形状および前記ダウンサンプリングされた重みテーブルから、計算処理装置を用いて補間によって再構築する段階も含む。
本発明の様々な実施形態によれば、非一時的なコンピュータ可読媒体が用意される。前記媒体は、プロセッサ上で実行されると、コンピュータ断層撮影装置に関連付けられた透視幾何学的形状データおよび重みテーブルをダウンサンプリングしかつ復元する命令を格納する。前記媒体は、前記コンピュータ断層撮影装置に関連付けられた透視幾何学的形状データを生成するための1つまたは複数の命令を格納する。前記媒体は、前記コンピュータ断層撮影装置に関連付けられた重みテーブルを生成するための1つまたは複数の命令をさらに格納する。格納するデータの量を減少させるために、前記重みテーブルおよび前記透視幾何学的形状データがダウンサンプリングされる。前記コンピュータ断層撮影装置に関連付けられた前記ダウンサンプリングされた透視幾何学的形状データおよび前記ダウンサンプリングされた重みテーブルが、記憶装置に格納される。前記媒体は、前記コンピュータ断層撮影装置を用いて走査された物体の断層撮影体積の再構築工程の間に、前記透視幾何学的形状データおよび前記重みテーブルを、それぞれ前記ダウンサンプリングされた透視幾何学的形状および前記ダウンサンプリングされた重みテーブルから、グラフィック処理ユニットを用いて補間によって再構築するための1つまたは複数の命令をさらに格納し、前記ダウンサンプリングされた重みテーブルが、前記GPUによってテクスチャ表面として格納される。
本発明の様々な実施形態によれば、1つのシステムが用意される。前記システムは、コンピュータ断層撮影装置に関連付けられた透視幾何学的形状データを生成するためのプロセッサを含む。前記プロセッサは、前記コンピュータ断層撮影装置に関連付けられた重みテーブルを生成するためにも使用される。格納するデータの量を減少させるために、前記重みテーブルおよび前記透視幾何学的形状データがダウンサンプリングされる。前記コンピュータ断層撮影装置に関連付けられた前記ダウンサンプリングされた透視幾何学的形状データおよび前記ダウンサンプリングされた重みテーブルが、記憶装置に格納される。前記プロセッサは、前記コンピュータ断層撮影装置を用いて走査された物体の断層撮影体積の再構築工程の間に、前記透視幾何学的形状データおよび前記重みテーブルを、それぞれ前記ダウンサンプリングされた透視幾何学的形状および前記ダウンサンプリングされた重みテーブルから、計算処理装置を用いて補間によって再構築するためにも使用される。
添付の図面は、本明細書の一部として組み込まれかつその一部をなすものであり、本明細書に記載された1つ又は複数の実施形態を図示し、その詳細な説明とともに実施形態を説明するものである。図面は次の通り。
(A) コンピュータ断層撮影を行うための代表的なCTスキャニングシステム100を示す。 (B) 図1Aのスキャニングシステム100に結合された計算装置を示す。 本発明の様々な実施形態による、重み計数を用いたコンピュータ断層撮影に関する三次元体積再構築を行うためのステップのフローチャートを示す。 (A) 重み値が、0度の投影角に関し再構築された体積のスライスにおいてどのように変化するかを示す。 (B) 重み値が、22.5度の投影角に関し再構築された体積のスライスにおいてどのように変化するかを示す。 (C) 重み値が、45度の投影角に関し再構築された体積のスライスにおいてどのように変化するかを示す。 (D) 重み値が、90度の投影角に関し再構築された体積のスライスにおいてどのように変化するかを示す。 (A) 本発明の様々な実施形態による重みテーブルを示す。 (B) 図4Aの重みテーブルに対応するダウンサンプリングされた重みテーブルを示す。
本発明は、逆投影を中央処理装置およびグラフィックス・アクセラレータで使用することにより、CT再構築の計算性能を向上させるための方法およびシステムを提供する。CTシステムは、選択されたスキャナの幾何学的形状を用いて、各反復投影の再構成体積の各平面におけるすべてのボクセルに関し、スキャナに関連付けられた透視幾何学的値および重み値を予め計算し、かつ、これらの結果を一連のテーブルに記憶する。本発明の様々な側面によれば、記憶されかつ逆投影時にロードされるテーブルのサイズを縮小でき、従って、メモリストレージ要件およびメモリアクセス帯域幅が減少する。本発明は、所与のスキャナに関連付けられた事前計算済みの透視幾何学的形状および重みテーブルをダウンサンプルし、次に、当該スキャナを用いて走査した物体の逆投影体積の計算時に、これらのテーブルを補間するための新規な方法を提供する。事前計算済みの重みテーブルを使用して計算要件を減少させ、次に、重みテーブルを補間するための付加的な計算を含めることは直観に反しかつ新規と考えられている。
再構築の間に重み値および透視幾何学的形状を正確に推定する必要があるため、テーブルのサイズは任意に縮小できない。また、本発明は、特定のコンピュータ断層撮影装置に関連付けられた重みテーブルおよび透視幾何学的形状が、再構築されたスライス内で小さな部分行列に関して徐々に変化するという新規な考察に基づく。この考察によれば、これらテーブルを一組のダウンサンプルされたテーブルとして格納することが奨励され、従って、断層撮影体積の再構築過程の間に、メモリストレージ要件およびメモリアクセス帯域幅要件が軽減される。断層撮影体積の再構築の間に、幾何および重みテーブルはそれ自体が必要に応じてある種の補間により再構築される。この革新的技法の価値の鍵は、多くのコンピュータアーキテクチャに関して、これらテーブルを再構築するのに必要な労力は、それらを元々の完全な形式で記憶するのに必要な労力より顕著に小さくなるということである。
図1Aは、代表的なCTスキャニングシステム100を示す。スキャニングシステム100では、供給源102は、z軸と位置合わせされた物体104の左側に位置している。検出器106は、物体104の右側に位置している。放射線ビームは、供給源102から検出器106に向けて発射される。供給源102から投射され物体104に吸収されない陽子は、検出器106に向かって送出されかつ検出器106に衝突し、検出器106上で物体102の像108を形成する。CT装置を構成する供給源102および検出器106は、物体104の周りを回転する。所定の複数角度および/または位置において、物体102の付加的な二次元投影像108が検出器106上に形成される。
図1Bに示したように、CTスキャニングシステム100の供給源102および検出器106、すなわちCT装置は計算装置120を用いて制御できる。計算装置120は、データを格納するための例えば記憶装置などのメモリ124と、後に詳述するグラフィック処理ユニット(GPU)などのプロセッサ122とを含むことができる。
計算装置120は、物体104の複数の二次元投影像108および各ボクセルに関する重み係数を用いて逆投影法により物体の三次元断層撮影体積を再構築する。重み係数は、供給源から物体の幾何学的形状の任意ボクセルへの距離と、検出器から物体の幾何学的形状の同じボクセルへの距離との比である。各ボクセルの重み係数は、予め計算し、例えばテーブル126の形式でメモリ124に保存しておくことができる。
図2は、本発明の様々な実施形態によるコンピュータ断層撮影装置に関連付けられたダウンサンプリング、透視幾何学的形状データおよび重みテーブルを行うためのステップのフローチャート200を示す。まず、コンピュータ断層撮影装置の透視幾何学的形状データが生成される(ステップ202)。ステップ204では、コンピュータ断層撮影装置の重みテーブルが生成される。これら幾何学的形状データおよび重みテーブルは、所与のCT装置に関連付けられている。これら幾何学的形状データおよび重みテーブルは走査する物体の関数として変化しない。従って、これら幾何学的形状データおよび重みテーブルは、同じCT装置を使った様々な物体の走査に再利用できる。このCT装置に関連付けられた走査速度および物理的な幾何学的形状が一定であれば、これら関連付けられた幾何学的形状データおよび重みテーブルも一定であり続ける。
上述のステップ202および204は計算が複雑であり、従って、非常にコストがかかる。しかし、特定のCT装置が一定の構成で動作する限りは、すなわちこのCT装置に関連付けられた走査速度および物理的な幾何学的形状が固定されている限り、幾何学的形状データおよび重みテーブルは、そのCT装置に関して一度だけ計算すればよい。
所与のCT装置に関連付けられた重みテーブルおよび透視幾何学的形状データがいったん計算されれば、これら重みテーブルおよび透視幾何学的形状データは、格納されるデータの量を減らすことによりダウンサンプリングされる(ステップ206)。ダウンサンプリングされた透視幾何学的形状データおよびダウンサンプリングされた重みテーブルは、記憶装置に格納される(ステップ208)。
走査された物体の断層撮影体積の再構築の間に、透視幾何学的形状データおよび重みテーブルが、それぞれダウンサンプリングされた透視幾何学的形状およびダウンサンプリングされた重みテーブルから、計算処理装置を用いて補間によって再構築される(ステップ210)。ステップ202ないし208は、実際の断層撮影体積の再構築工程が開始される前に行われる。ステップ210は、CT装置を用いて走査された各物体に関して、断層撮影体積の再構築工程時に連続的に繰り返される。
重み係数を示す重みテーブルは、物体104のすべての可能な(x,y,z)位置に関して計算できる。図1Aで示したように、xyz座標系には、中心点に示された物体104の1つのボクセルが存在している。具体的には、x軸およびy軸はz軸に直交する平面を形成する。uvs座標系には、中心点に物体104の投影像108が存在している。u軸はz軸に平行であり、v軸はz軸に直交する。重みテーブルの項目は次のように計算される。
Figure 0006101209
ここで、図1Aに示したベータ角(β)110は、y軸に対する供給源102の回転角である。xおよびy座標は、再構築中のスライス内の物理的位置である。座標(x, y)を備えた位置および回転角ベータ(β)に関して得られる重み係数は次の通り。
Figure 0006101209
DS0 112は、図1Aに示したように、供給源102から検出器106のu軸への直交光線の距離である。
図3A、3B、3C、および3Dはそれぞれ、重み値が、再構築された体積のスライスにおいてどのように変化するかを示す。投影角0度、22.5度、45度、および90度に関する重み行列または表面の平滑度を示す例が、それぞれ図3Aないし3Dに示されている。図3A、3B、3C、および3Dのそれぞれにおいて、垂直軸は、1つの投影角の二次元テーブルに格納された値を示す。これらの値は、1.7に近いピーク値から0.5に近い低値まで変化する。このグラフの他の2つの軸は、再構築された体積を通るスライスにおける位置(x,y座標)を表す。代表的な実施形態では、CTシステムにより使用される各投影角に関して、別の重みテーブルを計算し、格納できる。高分解能CTシステムは、各回転において数百あるいは数千にも上る別々の投影角を使用できる。例示目的で、4つの代表的な角度を図3Aないし3Dに示した。図3Aないし3Dに示した各グラフ内で、特定の値は、グラフの1つのセクションから別のセクションで大きく変化する。図3Aないし3Dに示したグラフは、平坦でなく浅い湾曲部を備えている。しかし、これらグラフは局所的に「平ら」である。すなわち、その勾配は、グラフにおける少数の隣接点にわたって大きく変化してはいない。従って、「変化率」の第1の微分はゼロに近い。
再構築計算の間には、重みテーブルが補間され、検出器上での投影ボクセルの位置および適用すべき重み因子を求めるために透視幾何学的形状を生成する。
例えば、システムは、それぞれの重み因子に関し、単精度フロート(4バイト)を用いて、合計400の投影で、1024 x 1024のボクセルを備えたスライスを再構築できる。対称性を利用しない場合は、このシステムでは、1.6ギガバイトのサイズの重みテーブルが得られる。この重みテーブルのサイズを減少させれば、メモリストレージおよび帯域幅要件を顕著に減少できる。従って、システム性能およびコストにおいて顕著な向上が得られる。重みテーブルをxおよびy軸において2の係数でダウンサンプリングことにより、重みテーブルのメモリストレージおよびアクセス要件が4分の一に減少する。
図4Aは、代表的な重みテーブル400の小さな(5x5)部分を示す。この重みテーブルは、双一次、バイキュービック、または周波数領域におけるダウンサンプリングなどを含むがそれらに限定されない方法を用いてダウンサンプリングできる。図4Bは、図4Aに示した完全な重みテーブル400に対応した代表的なダウンサンプリングされた重みテーブル402を示す。代表的なダウンサンプリングされた重みテーブル402は、最近隣選択法(原語: nearest neighbor selection method)を用いて計算される。逆投影計算の実行時に、このダウンサンプリングされた重みテーブルは、双一次リサンプリング方式などの補間方式を用いて、完全な、すなわちダウンサンプリングされていない重みテーブルに再計算できる。例えば、この重みテーブルが2の係数でダウンサンプリングされていれば、双一次補間すなわち双一次リサンプリングの係数は、0および1または0.5および0.5である。すなわち、ボクセル(Y=2, X=2)の重みW22の推定値を生成するには、ダウンサンプリングされた行列402から次のように行列値を組み合わせる必要があり、すなわち、W22 ~= 0.5 * (0.5*W11 + 0.5*W13) + 0.5 *
(0.5*W31+0.5*W33)である。
上述の例は、例示目的のみで記載したものであり、限定的に解釈すべきではない。
ダウンサンプリングの量は、テーブルの度数内容(原語: frequency content)と再構築で必要とされる忠実度により制限される。このダウンサンプリングは、2の係数に限定されるものでなく、両方の次元において同一である必要もない。1024x1024ボクセルのスライスを生成する代表的なシステムでは、重み行列が、各次元において16でダウンサンプリングされれば、1024x1024の値の重みテーブルは64x64配列に縮小できる。こうしたダウンサンプリングは、重みテーブルの合計メモリ要件を256の係数で縮小する。
重みテーブルのローディングおよびその後の補間は様々な方法で実行可能であり、重みテーブルの異なる格納モデルの恩恵を受けることができる。本発明の実施形態によれば、逆投影法はグラフィック処理ユニット(GPU)で実装できる。
代表的な実施形態は、CUDAソフトウェア実装を用いたnVIDIA Fermiアーキテキクチャに基づいている。この例では、重み行列は、最近隣選択法を用いてxおよびy次元で2の係数でダウンサンプリングされる。逆投影時に更新されるスライスは1024x1024のサイズであり、従って格納されるダウンサンプリングされた重みテーブルは512x512のサイズである。スライス内のピクセル位置ix, iy (ixおよびiyは両方とも0から1023の範囲を備えている)を更新する際には、重みテーブルのインデックスはix/2およびiy/2である。従って、分数テーブルアドレスix/2, iy/2を用いて適用される重み因子は、重み因子Wt(ix,iy) = Wt_downsample(ix/2, iy/2)である。
この重みテーブルは、計算装置のグローバルデバイスメモリまたはテクスチャメモリに格納できる。テクスチャメモリは、自動補間を行うというGPUハードウェアの独自の機能を与える。すなわち、重みテーブルの位置は、整数アドレスでなく分数アドレスに選択できる。
GPUハードウェアは、最近隣または双一次補間を用いて非常に効果的に補間を実行できる。ダウンサンプリングされたテクスチャを補間して完全な解像度行列を満たす方法を示す代表的なコードを後に示す。nVidia社のCUDA Cでは、プログラマーは、呼び出されると通常のC関数のように一度だけでなくN個の異なるCUDAスレッドによって並列でN回実行されるカーネル(CUDAカーネル)と呼ばれるC関数を定義できる。下記のinterpolateWeight関数はそうしたCUDAカーネルである。変数downsampledWeightTextureはテクスチャ表面への参照であり、full_resoution_weightは生成される完全な解像度重み行列であり、down-sampleはダウンサンプル係数、sizeは完全な解像度重み行列の幅および高さである。blockDim、blockIdx、およびthreadIdxは、各CUDAスレッドにより実行されるCUDAカーネルが、どのデータを処理するかを決定できるようにする特別な構造体である。位置uおよびvは浮動小数点値であり、従って、分数すなわち非整数位置においてテクスチャ表面として格納されたテーブルをサンプリングする。
// 2D浮動小数点テクスチャは二次元テクスチャを定義。.
texture<float, 2, cudaReadModeElementType> downsampledWeightTexture;
// Kernel definition
__global__ void interpolateWeight(float* full_resolution_weight, int downsample, int size)
{
// 正規化テクスチャ座標を計算。
unsigned int x = blockIdx.x * blockDim.x + threadIdx.x;
unsigned int y = blockIdx.y * blockDim.y + threadIdx.y;
float u = x / (float)downsample;
float v = y / (float)downsample;

if(y<size && x<size)

full_resolution_weight[y * size + x] = tex2D(downsampledWeightTexture, u,
v);
}
}
GPUにおける実施形態はテクスチャ表面に限定されない。本発明の様々な実施形態によれば、各処理カーネルがダウンサンプリングされたテーブルの四隅をロードし、CUDAカーネル内で補間を計算することが有利となることがある。
汎用プロセッサを使って同様のアプローチをとることができる。その場合、テーブルのロードは4タプルとして行うことができ、ベクトルプロセッサ(すなわち、AltivecまたはIntel SSE)を用いて補間値を計算できる。
どの実施形態を使用した場合も、同一の根本的な革新から利益がもたらされるが、その革新とは、重みテーブルのダウンサンプリングされた表現の格納ならびに逆投影時のこれらテーブルを補間することによる帯域幅およびメモリストレージ量の減少である。
重み因子およびボクセルの検出器への投影を定義する透視幾何学的形状が、後述する予め計算された配列により定義される再構築を考慮する。テーブルは一組のz平面からなんらかの周期性があると想定されるが、この例に関してはその周期性は示されていない。
Xi 0 <= i < Nx は、所与のスライスに関するボクセルのXインデックスを定義し、ここで、Nxはxにおけるボクセル数である。
Y j 0 <= j < Ny は、所与のスライスに関するボクセルのYインデックスを定義し、ここで、Nyはyにおけるボクセル数である。
Zk 0 <= k < Nz は、Zインデックス(またはスライス)を定義し、ここで、Nzはzにおけるボクセル(またはスライス)数である。
V l 0<= l< Nvは、検出器における行インデックスを定義し、ここで、Nvは検出器の行の数である。
U m 0<= m < Nuは、検出器における列インデックスを定義し、ここで、Nuは検出器の列の数である。
αn 0<= n < Nαは、投影インデックスを定義し、ここで、Nαは1つのスライスに寄与する投影の数である。
Vol(Zk, Yj, Xi)は再構築される体積である。
P(z, α,Um, Vl)は、u,v グリッドにおける検出器からの二次元生データである。
U(z, α, y, x)は、各エレメントが、U(Zk, αn, Y j, Xi)を選択することによって得られる位置を定義する。所与のZk, αn,に関し、この位置は二次元テーブルである。
V(Zk, αn, Y j, Xi)を選択することによって各エレメントが得られる位置V(z, α, y,
x)。所与のZk, αn,に関し、この位置は二次元テーブルである。
W(Zk, αn, Y j, Xi)を選択することによって各エレメントが得られる位置W(z, α, y,
x)。所与のZk, αn,に関し、この位置は二次元テーブルである。
上述のように、所与のZ平面および所与の投影αnに関し、二次元テーブルを格納できる。処理の間に補間できるのはこれら二次元テーブルである。
逆投影ステップの実行は次のステップに限定される。
For each (Zk)
For
each project(αn)
For each Yj

For each Xi

Vol(Zk, Yj, Xi) = Vol(Zk, Yj, Zi)
+ W(Zk, αn, Yj, Xi) * P(Zk, αn, V(z, α,
y, x), U(z, α, y, x));

End

End
End
End
この例では、重み行列ならびに検出器上の投射の位置UおよびVは、完全な解像度行列に格納される。しかし、本発明は、処理の前にこれらテーブルをダウンサンプリングし、逆投影の内側ループの間に補間することを提案する。従って、例えば、ボクセルの検出器上への投影位置および重みは、すべて4でダウンサンプリングされていると想定する。ダウンサンプリングされた行列をW4、V4、およびU4で表現する。
すると、処理ステップは次のようになる。
For each (Zk)
For
each project(αn)
// 二次元位置および重み行列への2つの参照(ポインタ)を得る。
WW4 = W4(Zk, αn, :, :); # 二次元行列の開始点のみを得る。
UU4 = U4(Zk, αn, :,: )
VV4 = V4(Zk, αn, :,: )
For each Yj

For each Xi

U = interp(WW4, Yj, Xi);

V = interp(UU4, Yj, Xi);

WT = interp(WW4, Yj,Xi);

Vol(Zk, Yj, Xi) = Vol(Zk, Yj, Zi)
+ WT * P(Zk,
αn, V, U);

End

End
End
End
補間の詳細は、アーキテクチャと格納方法に依存する。
実施形態および実例の上記説明は例示および解説を意図したものであり、すべてを網羅することを意図してはいない。修正および変更は上述の教示を参照すれば可能であるし、本発明を実行することで発見されることもあろう。従って、本発明は、次に添付された特許請求の範囲内に収まるすべての実施形態および均等物を含むことが意図されている。
本明細書に記載された1つまたは複数の実施形態および実例は、多数の異なる形式のソフトウェアおよびハードウェアで実装できることは明らかなはずである。本明細書に記載された実施形態を実装するために使用されるソフトウェアコードおよび/または専用ハードウェアは、本発明を限定するものではない。従って、実施形態の動作および挙動は、これら特定のソフトウェアコードおよび/または専用ハードウェアを参照することなく記載されている。すなわち、本明細書の記載に基づけばこれら実施形態を実装するためのソフトウェアおよび/またはハードウェアを設計できることは理解されるはずである。

Claims (16)

  1. コンピュータ断層撮影装置に関連付けられた透視幾何学的形状データおよび重みテーブルをダウンサンプリングしかつ復元するためのコンピュータ実装方法であって、
    前記コンピュータ断層撮影装置に関連付けられた透視幾何学的形状データを生成する段階と、
    前記コンピュータ断層撮影装置に関連付けられた重みテーブルを生成する段階と、
    格納するデータの量を減少させるために、前記重みテーブルおよび前記透視幾何学的形状データをダウンサンプリングする段階と、
    前記コンピュータ断層撮影装置に関連付けられた前記ダウンサンプリングされた透視幾何学的形状データおよび前記ダウンサンプリングされた重みテーブルを記憶装置に格納する段階と、
    前記コンピュータ断層撮影装置を用いて走査された物体の断層撮影体積の再構築工程の間に、前記透視幾何学的形状データおよび前記重みテーブルを、それぞれ前記ダウンサンプリングされた透視幾何学的形状および前記ダウンサンプリングされた重みテーブルから、計算処理装置を用いて補間によって再構築する段階とを含み、
    前記計算処理装置がグラフィック処理ユニット(GPU)であり、
    前記ダウンサンプリングされた重みテーブルが、前記GPUによってテクスチャ表面として格納される、コンピュータ実装方法。
  2. 前記断層撮影体積の再構築工程の間に、前記透視幾何学的形状データおよび前記重みテーブルを用いて前記物体の三次元断層撮影体積を生成する段階をさらに含む、請求項1に記載の方法。
  3. 前記三次元断層撮影体積は、逆投影法、代数的再構築手法(ART)、または最尤推定-期待値最大化(MLEM)法を用いて生成される、請求項2に記載の方法。
  4. 前記重みテーブルは二次元テーブルである、請求項1に記載の方法。
  5. 前記重みテーブルは、コンピュータ断層撮影のために走査される前記物体の幾何学的形状内部の各ボクセルに関する重みを含む、請求項1に記載の方法。
  6. 任意のボクセルに関する前記重みは、前記ボクセルから前記コンピュータ断層撮影に用いられるエネルギー源への距離と、前記ボクセルから前記コンピュータ断層撮影に用いられる検出器への距離との比に等しい、請求項5に記載の方法。
  7. 前記ダウンサンプリングされた重みテーブルは、双一次ダウンサンプリング、バイキュービックダウンサンプリング、または周波数領域におけるダウンサンプリングの1つまたは複数を用いて前記重みテーブルから計算される、請求項1に記載の方法。
  8. プロセッサであって、
    コンピュータ断層撮影装置に関連付けられた透視幾何学的形状を生成し、
    前記コンピュータ断層撮影装置に関連付けられた重みテーブルを生成し、
    格納するデータの量を減少させるために、前記重みテーブルおよび前記透視幾何学的形状データをダウンサンプリングし、
    前記コンピュータ断層撮影装置に関連付けられた前記ダウンサンプリングされた透視幾何学的形状データおよび前記ダウンサンプリングされた重みテーブルを記憶装置に格納し、
    前記コンピュータ断層撮影装置を用いて走査された物体の断層撮影体積の再構築工程の間に、前記透視幾何学的形状データおよび前記重みテーブルを、それぞれ前記ダウンサンプリングされた透視幾何学的形状および前記ダウンサンプリングされた重みテーブルから補間によって再構築するプロセッサを含み、
    前記プロセッサがグラフィック処理ユニット(GPU)であり、
    前記ダウンサンプリングされた重みテーブルが、前記GPUによってテクスチャ表面として格納される、システム。
  9. 前記プロセッサが、前記断層撮影体積の再構築工程の間に、前記透視幾何学的形状データおよび前記重みテーブルを用いて前記物体の三次元断層撮影体積をさらに生成する、請求項8に記載のシステム。
  10. 前記三次元断層撮影体積は、逆投影法、代数的再構築手法(ART)、または最尤推定-期待値最大化(MLEM)法を用いて生成される、請求項9に記載のシステム。
  11. 前記重みテーブルは二次元テーブルである、請求項8に記載のシステム。
  12. 前記重みテーブルは、コンピュータ断層撮影のために走査される前記物体の幾何学的形状内部の各ボクセルに関する重みを含む、請求項8に記載のシステム。
  13. 任意のボクセルに関する前記重みは、前記ボクセルから前記コンピュータ断層撮影に用いられるエネルギー源への距離と、前記ボクセルから前記コンピュータ断層撮影に用いられる検出器への距離との比に等しい、請求項12に記載のシステム。
  14. 前記ダウンサンプリングされた重みテーブルは、双一次ダウンサンプリング、バイキュービックダウンサンプリング、または周波数領域におけるダウンサンプリングの1つまたは複数を用いて前記重みテーブルから計算される、請求項8に記載のシステム。
  15. コンピュータ上で実行されると、コンピュータ断層撮影装置に関連付けられた透視幾何学的形状データおよび重みテーブルをダウンサンプリングしかつ復元する命令を格納するための非一時的なコンピュータ可読媒体であって、
    前記コンピュータ断層撮影装置に関連付けられた透視幾何学的形状データを生成し、
    前記コンピュータ断層撮影装置に関連付けられた重みテーブルを生成し、
    格納するデータの量を減少させるために、前記重みテーブルおよび前記透視幾何学的形状データをダウンサンプリングし、
    前記コンピュータ断層撮影装置に関連付けられた前記ダウンサンプリングされた透視幾何学的形状データおよび前記ダウンサンプリングされた重みテーブルを記憶装置に格納し、
    前記コンピュータ断層撮影装置を用いて走査された物体の断層撮影体積の再構築工程の間に、前記透視幾何学的形状データおよび前記重みテーブルを、それぞれ前記ダウンサンプリングされた透視幾何学的形状および前記ダウンサンプリングされた重みテーブルから、グラフィック処理ユニット(GPU)を用いて補間によって再構築するための1つまたは複数の命令を格納し、
    前記ダウンサンプリングされた重みテーブルが、前記GPUによってテクスチャ表面として格納される、非一時的なコンピュータ可読媒体。
  16. 前記断層撮影体積の再構築工程の間に、前記透視幾何学的形状データおよび前記重みテーブルを用いて前記物体の三次元断層撮影体積を生成するための1つまたは複数の命令をさらに格納する、請求項15に記載の非一時的なコンピュータ可読媒体。
JP2013550621A 2011-01-21 2012-01-20 コンピュータ断層撮影(ct)のための逆投影の最適化実装 Active JP6101209B2 (ja)

Applications Claiming Priority (5)

Application Number Priority Date Filing Date Title
US201161435107P 2011-01-21 2011-01-21
US61/435,107 2011-01-21
US13/354,587 2012-01-20
PCT/US2012/022051 WO2012100177A1 (en) 2011-01-21 2012-01-20 Optimized implementation of back projection for computed tomography (ct)
US13/354,587 US8494213B2 (en) 2011-01-21 2012-01-20 Optimized implementation of back projection for computed tomography (CT)

Publications (2)

Publication Number Publication Date
JP2014506491A JP2014506491A (ja) 2014-03-17
JP6101209B2 true JP6101209B2 (ja) 2017-03-22

Family

ID=46516118

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2013550621A Active JP6101209B2 (ja) 2011-01-21 2012-01-20 コンピュータ断層撮影(ct)のための逆投影の最適化実装

Country Status (6)

Country Link
US (1) US8494213B2 (ja)
EP (1) EP2666122A4 (ja)
JP (1) JP6101209B2 (ja)
AU (1) AU2012207195B2 (ja)
CA (1) CA2825704A1 (ja)
WO (1) WO2012100177A1 (ja)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6280851B2 (ja) * 2014-09-30 2018-02-14 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 放射線断層撮影装置及びプログラム
CN104574509A (zh) * 2015-01-26 2015-04-29 上海交通大学 一种由投影重建物体三维图像的方法
US9626779B1 (en) * 2015-10-20 2017-04-18 International Business Machines Corporation Efficient back-projection operation using precomputed table
WO2019096600A1 (en) * 2017-11-14 2019-05-23 Koninklijke Philips N.V. Single ct backprojector with one geometry calculation per voxel for multiple different types of projection data
CN113017662B (zh) * 2021-01-28 2022-06-14 明峰医疗系统股份有限公司 一种ct图像的混叠伪影去除方法及系统、ct扫描仪

Family Cites Families (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4135247A (en) * 1977-08-15 1979-01-16 Siemens Aktiengesellschaft Tomography signal processing system
US4809065A (en) * 1986-12-01 1989-02-28 Kabushiki Kaisha Toshiba Interactive system and related method for displaying data to produce a three-dimensional image of an object
US5982845A (en) * 1998-04-07 1999-11-09 Picker International, Inc. Forward projection and backprojection processor
EP1700270B1 (en) * 2003-12-16 2010-06-23 Philips Intellectual Property & Standards GmbH Imaging method with filtered back projection
DE102004022332A1 (de) * 2004-05-06 2005-12-08 Siemens Ag Verfahren zur post-rekonstruktiven Korrektur von Aufnahmen eines Computer-Tomographen
US7376255B2 (en) 2004-06-23 2008-05-20 General Electric Company System and method for image reconstruction
US7355598B2 (en) * 2004-06-25 2008-04-08 Siemens Medical Solutions Usa, Inc. System and method for fast generation of high-quality maximum/minimum intensity projections
US7149276B2 (en) * 2004-07-14 2006-12-12 Kabushiki Kaisha Toshiba System, method, and computer program product that corrects measured data
US7778392B1 (en) * 2004-11-02 2010-08-17 Pme Ip Australia Pty Ltd Method of reconstructing computed tomography (CT) volumes suitable for execution on commodity central processing units (CPUs) and graphics processors, and apparatus operating in accord with those methods (rotational X-ray on GPUs)
JP5269283B2 (ja) * 2005-04-06 2013-08-21 株式会社東芝 医用画像診断装置及び画像再構成方法
CN100444798C (zh) * 2005-07-22 2008-12-24 清华大学 一种用于ct重建的旋转对称体素离散化方法
CN101416072B (zh) * 2005-09-02 2012-03-14 皇家飞利浦电子股份有限公司 用于计算机断层造影成像的改进的重组
WO2007034342A2 (en) * 2005-09-26 2007-03-29 Koninklijke Philips Electronics, N.V. Iterative reconstruction with enhanced noise control filtering
JP4414420B2 (ja) * 2006-10-27 2010-02-10 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー X線断層撮影装置およびアーチファクトの低減方法
US20090195541A1 (en) 2008-02-05 2009-08-06 Rambus Inc. Rendering dynamic objects using geometry level-of-detail in a graphics processing unit
US7852977B2 (en) * 2008-09-11 2010-12-14 Samplify Systems, Inc. Adaptive compression of computed tomography projection data
US7829856B2 (en) * 2009-03-31 2010-11-09 General Electric Company Apparatus and methods for determining a system matrix for pinhole collimator imaging systems
DE112009005019B4 (de) * 2009-06-30 2022-02-03 Analogic Corp. Effizienter quasi-exakter 3D Bildrekonstruktionsalgorithmus für CTScanner

Also Published As

Publication number Publication date
EP2666122A4 (en) 2018-01-24
US20120189158A1 (en) 2012-07-26
EP2666122A1 (en) 2013-11-27
CA2825704A1 (en) 2012-07-26
AU2012207195B2 (en) 2016-09-15
US8494213B2 (en) 2013-07-23
JP2014506491A (ja) 2014-03-17
AU2012207195A1 (en) 2013-09-05
WO2012100177A1 (en) 2012-07-26

Similar Documents

Publication Publication Date Title
Pratx et al. Fast, accurate and shift-varying line projections for iterative reconstruction using the GPU
Rahmim et al. Statistical dynamic image reconstruction in state-of-the-art high-resolution PET
US9801591B2 (en) Fast iterative algorithm for superresolving computed tomography with missing data
US9002090B2 (en) Forward projection apparatus
JP6081477B2 (ja) Ct装置及びct画像生成方法
JP2008006288A (ja) 繰り返し式画像再構成のシステム及び方法
JP6101209B2 (ja) コンピュータ断層撮影(ct)のための逆投影の最適化実装
JPH0910202A (ja) 診断撮像方法と装置およびヘリカル部分コーンビームデータからの画像再構成
Liu et al. GPU-based branchless distance-driven projection and backprojection
Keck et al. GPU-accelerated SART reconstruction using the CUDA programming environment
Kim et al. Fully 3D iterative scatter-corrected OSEM for HRRT PET using a GPU
Kole et al. Evaluation of accelerated iterative x-ray CT image reconstruction using floating point graphics hardware
Sunnegårdh et al. Regularized iterative weighted filtered backprojection for helical cone‐beam CT
US9495770B2 (en) Practical model based CT construction
US7272205B2 (en) Methods, apparatus, and software to facilitate computing the elements of a forward projection matrix
Yendiki et al. A comparison of rotation-and blob-based system models for 3D SPECT with depth-dependent detector response
Johnston et al. GPU-based iterative reconstruction with total variation minimization for micro-CT
Vintache et al. Iterative reconstruction for transmission tomography on GPU using Nvidia CUDA
US12045916B2 (en) Stochastic backprojection for 3D image reconstruction
Steckmann et al. High performance cone-beam spiral backprojection with voxel-specific weighting
Chen et al. Accelerating ring artifact correction for flat-detector CT using the CUDA framework
Stsepankou et al. Real-time 3D cone beam reconstruction
Huang et al. Rapid prejudgment of reconstructed object volume and its adaptive reconstruction for industrial cone-beam CT
Gao et al. Distributed computation of linear inverse problems with application to computed tomography
Azhari et al. Basic Principles of Tomographic Reconstruction

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20141212

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20150831

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20151006

A02 Decision of refusal

Free format text: JAPANESE INTERMEDIATE CODE: A02

Effective date: 20160706

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20161027

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A821

Effective date: 20161121

A911 Transfer to examiner for re-examination before appeal (zenchi)

Free format text: JAPANESE INTERMEDIATE CODE: A911

Effective date: 20161228

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: 20170131

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20170224

R150 Certificate of patent or registration of utility model

Ref document number: 6101209

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250