JP2014505864A - 動的な多次元画像データの再構成 - Google Patents

Abstract

一連の再構成時点で動的なサンプルを表す多次元データセットを再構成する方法が開示される。多次元データセットは、静的な要素と動的な要素とを有する。この方法は、動的なサンプルの複数の投影画像を取得し、取得された投影画像から多次元データセットの静的な要素を再構成し、動的なサンプルの更なる複数の投影画像を取得し、サンプルについての事前情報を使用して、更なる複数の取得された投影画像から各再構成時点での多次元データセットの動的な要素を再構成することを有する。多次元データセットは、各再構成時点での静的な要素及び動的な要素の和である。

Description

本発明は、概して断層撮影画像化に関し、特に、断層撮影スキャンからの動的な多次元画像データの再構成に関する。

X線コンピュータ断層撮影法（CT：computed tomography）は、単一の“スキャン（scan）”を構成する異なる視野角の範囲から多次元サンプル（標本）の複数の１次元又は２次元（1D又は2D）放射線写真を取得することにより実行される。1D又は2D投影画像のこのセットから、サンプル内のX線の線減衰係数の2次元又は3次元空間分布を示すサンプルの多次元（2D又は3D）画像が再構成され得る。高分解能のX線CT（マイクロCT）を使用した複雑な内部構造の非破壊検査は、材料科学、生物学、地質学及び石油工学のような分野で急速に標準的な技術になっている。歴史的に、CTは、単に静的なサンプルを画像化するために使用されてきている。例外は医療用の画像化であり、特に心臓又は肺のような動く臓器を画像化するときに、患者の動きがしばしば回避できない。医療の状況でのダイナミックCTは、心拍のような動きを訂正し、時間変化する要素を除去することにより、高品質な静的な画像を形成することを試みる画像化技術を示す。

これに対して、ほとんどのマイクロCT画像化では、3Dサンプルの動的なふるまい（時間変化）が真に関心のあるものであり、必ずしも意図しない動き又は周期的な動きの結果であるとは限らない。例えば、多孔質材料の内部の１つの不混和性の他のものによる流体の置換は、基礎となる物理の複雑性のため、また、標準的な実験が微小規模のプロセスについて非常に少ないことしか明らかにしていないため、地質学において特に困難な問題である。多相の置換は、石油生産の中核である。この理由は、地層において水が油に置換する方法は、油がその地層から経済的に抽出できるか否か及びどのように抽出できるかを決定するからである。同相の４次元（4D）実験データ（時間軸上の3D）は、取得するのに極めて高コストであり、失望するほど少ない情報しか戻さない。モデル化の研究が安価であり、多くの見識を提供するが、真の予測力を欠いている。実験とモデルとの間の微小規模での比較は、モデル化が有用であるためには非常に必要になる。ダイナミックマイクロCTは、研究室の条件でのこのような4D実験データを取得するための基本的に適切な様式である。

放射線画像セットでCT再構成を実行する従来の方法は、フィルターバック投影（FBP：filtered backprojection）と、フーリエ反転と、様々な反復方式（代数的再構成法（ART：algebraic reconstruction technique）、同時反復再構成法（SIRT：simultaneous iterative reconstruction technique）及び関係する同時代数的再構成法（SART：simultaneous algebraic reconstruction technique）等）とを含む。このような技術の全ては、(i)サンプルが静的であり、(ii)サンプル内の関心のある構造が完全に各放射線写真の視野内にあることを仮定している。サンプルが取得中に変化した場合、放射線写真は、互いに一貫性がなくなり、再構成された画像のアーチファクト（artifact)及び／又はぶれをもたらす。実際に、このことは、従来のCT画像化が、サンプルが放射線写真の全セットを取得する時間の間に事実上静的である状況に制限されることを意味する。

CT再構成の問題は、高周波数の実験雑音に関して少し不安定であり、N³グリッドのボリュームエレメント（ボクセル）上で3D画像を正確に再構成するために、約π/2の視野角での放射線写真が必要になることが証明されている。各放射線写真の取得時間はN²に比例するため、従来のCT再構成を使用した最大の実現可能な時間分解能はN³に比例する。言い換えると、従来のCTを使用して、サンプルが基本的に静的に留まらなければならない時間量は、所望の空間分解能に比例して増加する。

研究室に基づくCTシステムでは、空間分解能の増加は、小さい領域からX線を放射する源を使用することを意味し、電子ビームがターゲットの物質の小さい領域に集中されなければならないことを意味する。このことは、ビーム出力を基本的に制限する。この理由は、小さすぎる領域に集中する大きすぎるエネルギーはターゲットの物質を気化させるからである。また、このことは、許容可能な信号対雑音比（SNR：signal-to-noise ratio）で単一の放射線写真を取得するために必要な時間量に低い制限を課す。従って、高分解能の研究室に基づくCTスキャンは、典型的には4〜15時間を要し、現在関心のある動的に変化するサンプルを画像化するには許容できない時間分解能である。

既存の構成の１つ以上の欠点を実質的に克服すること又は少なくとも改善することが、本発明の目的である。

従来の研究室に基づくCT画像化システムを使用した動的な（時間変化する）2D又は3Dサンプルのダイナミック（2D+時間又は3D＋時間）CT画像化の方法が開示される。開示の方法は、小さい数の取得された投影画像から動的な画像の安定した再構成を可能にし、これにより従来の方法に対して時間分解能を改善するために、画像化されているサンプルについての事前情報を利用する。

本発明の第１の態様によれば、一連の再構成時点で静的な要素と動的な要素とを有する動的なサンプルを表す多次元データセットを再構成する方法が提供される。この方法は、動的なサンプルの複数の投影画像を取得し、取得された投影画像から多次元データセットの静的な要素を再構成し、動的なサンプルの更なる複数の投影画像を取得し、動的なサンプルについての事前情報を使用して、更なる複数の投影画像から各再構成時点での多次元データセットの動的な要素を再構成することを有する。多次元データセットは、各再構成時点での静的な要素及び動的な要素の和である。

本発明の第２の態様によれば、複数の取得時点及び対応する視野角で取得されたサンプルの投影画像のセットから一連の再構成時点での動的なサンプルを表す一連の多次元画像を再構成する方法が提供される。この方法は、各再構成時点で、視野角における再構成時点での多次元画像の現在の推定を投影し、複数の投影を形成し、投影及び投影画像の系列から差分画像を形成し、ただし、系列は、再構成時点の周辺の取得時点で取得された連続する投影画像を有し、サンプルを通じて投影経路長により差分画像を正規化し、正規化された差分画像を逆投影し、再構成時点での多次元画像の現在の推定に逆投影を加算することを有する。

本発明の他の態様も開示される。

本発明の実施例が記載される座標系での典型的な円錐ビームCT画像化幾何学配置を示す図 実施例に従って図１の画像化幾何学配置を使用して取得された画像のセットから動的に変化するサンプルの線減衰係数の動的な要素を再構成する方法を示すフローチャート 図２の方法が実施され得る汎用コンピュータシステムの概略ブロック図 図２の方法が実施され得る汎用コンピュータシステムの概略ブロック図 実験データセットで図２の方法を実行した結果を示す図

本発明の少なくとも１つの実施例について図面を参照して説明する。

１つ以上の添付図面において同じ参照符号を有するステップ及び／又は機能を参照した場合、これらのステップ及び／又は機能は、逆の意図が示されない限り、この説明の目的で同じ機能又は動作を有する。

ダイナミックCT再構成の開示の方法は、サンプルの静的な特徴を“取り除き”、取得時点と次の時点との間に生じる比較的小さい変化に集中する。直感的に、これは、動的な信号ができるだけ小さいデータを使用して符号化される動画像符号化に類似する。再構成の問題におけるデータ圧縮方式（“圧縮センシング（compressed sensing）”と呼ばれる）の使用は、音響数学の基礎（sound mathematical foundation）を有する。圧縮センシング（CS）法は、最小の空間的変化を示す静的なサンプルのCT再構成のために開発された。このような“CS-CT”法は、コスト関数が(i)解と測定データとの間のずれと(ii)解の“全変動（total variation）”（すなわち、グラディエントのLIノルム）との適切に重み付けられた組み合わせである最適化問題としてCT再構成を扱う。このコスト関数は、サンプルについての事前情報と一致する全ての画像で最小化される。

離散断層撮影法（DT：discrete tomography）では、静的なサンプルが数個の（典型的には２つの）濃度値のみを使用して表され得ることを事前に想定する。DT法は、別法では必要とされるものよりもかなり少ない放射線写真から再構成を可能にし、スキャン時間及び／又はX線照射線量の比例した低減をもたらす。開示の方法はまた、動的なサンプルについての適切な事前情報を利用する。これにより、再構成の問題は、スキャン時間と再構成の空間分解能との間の比例関係を壊すように変更される。これにより、約１オーダーの大きさでダイナミックマイクロCT画像化の改善した時間分解能を可能にする。

開示の方法は、事前情報から導かれた適切な制限を用いて、最適化問題として再構成の問題を定式化する。従って、開示の方法は、DTとCSとの双方の要素を組み込む。開示の方法は、(i)複雑であるが静的な特徴と(ii)事前情報が定式化され得る動的な特徴とで構成されたサンプルに適用可能である。

図１は、本発明の実施例が定式化される座標系での典型的な円錐ビームCT画像化幾何学配置100を示している。円錐ビームCT画像化幾何学配置100は、従来の研究室に基づくCT 3D画像化システムをモデル化するのに適している。2D画像化幾何学配置は、2Dサンプル（r₃=一定）と1（以上）のピクセルの高さの検出器とを用いた3D画像化幾何学配置100の特別な場合として定式化されてもよい。記載の再構成法は、平面ビーム、ファンビーム（fan-beam）、らせん（helical）等の十分にうまく動作する投影及び逆投影演算子を用いた他のCT画像化幾何学に適用されてもよい。

3Dデカルトサンプル空間座標系r=(r₁,r₂,r₃)は、画像化されている2D又は3Dサンプル110に対して固定されており、中心にある。サンプル110は、その線減衰係数μ(r,t)を使用してモデル化される。ただしtは時間である。サンプル110は、最初に（すなわち、t<0では）静止状態にある。t=0において、動的なプロセスが開始される（例えば、ポンプがオンになる）。概念的に言えば、開示の方法は、時間毎に生じる変化を分離し、これらの変化を再構成する。１つの実装では、サンプルはμ_s(r)としてモデル化された静的な要素と、μ_d(r,t)としてモデル化された絶えず変化する動的な要素とに分離される。

１つの実装では、線減衰係数μ(r,t)についての事前情報は３つの仮定を有する。

1.動的な要素μ_d(r,t)が、（例えば、非圧縮性流体、二相流、反応流又は岩石（rock）の圧縮（compression）等で）数個の濃度値のみを使用して正確に表され得る。微孔性岩石の非反応性流体の流れの場合、濃度値の数は、流れの存在／不在を表す２つである。N相反応流の場合、岩石の浸食を表すために更なる負の濃度値が追加されてもよい。この新たな濃度値は、以前に岩石で満たされた領域にのみ存在してもよい。岩石圧縮の用途では、仮定1は、１つの再構成時点から次の時点への濃度値及び質量の保存を実現するように定式化される。

2.動的な要素μ_d(r,t)の瞬間的な変化が小さい。微孔性岩石を通じた非反応性流体又は反応性流体の流れの場合、瞬間的な変化は空間的に局所化される。他の用途（例えば、岩石圧縮では、瞬間的な変化はサンプルを通じて生じる。このような用途では、仮定2の適切な定式化は、サンプルの特徴が遠くに動かないことである。

3.動的な要素μ_d(r,t)の支持情報（例えば、細孔又はサンプル境界）が、静的な要素μ_s(r)から導かれ得る。不浸透性の微孔性の静的な骨組（scaffold）を通じた非反応性流体の流れの場合、この仮定は、骨組の細孔の空間（pore-space）（静的な要素μ_s(r)の支持領域の要素）が流体（動的な要素μ_d(r,t)）の支持領域であるように定式化される。反応流の用途では、動的な要素μ_d(r,t)の支持領域は、骨組が浸食されると共に時間と共に増加する。この仮定は、流体の流れの既知の特性に整合したレートで動的な要素μ_d(r,t)の支持領域をゆっくり拡大することにより定式化されてもよい。岩石圧縮の用途では、サンプルの構造、質量及び境界は、静的な要素μ_s(r)から導かれてもよい。

サンプル110は、サンプルの座標系rにおける位置s(θ)=(R₁cosθ,R₁sinθ,0)で、サンプルからの距離R₁においてI_inの強度（サンプルへの入射）の微小焦点X線源から発散する放射で照射される。θは、水平面における源とサンプルとの間の“視野角”である。（シンクロトロンで典型的に存在する平行平面のスキャン幾何学配置は、R₁が無限大に向かう円錐ビーム幾何学配置100の限定的な場合である。）
X線強度は、サンプル110の“下流”のR₂の距離において検出面の2D位置検知形検出器130により測定される。検出器130は、源120に対して固定されている。x=(x1,x2)は、源120から検出器130への軸に垂直な起点r=0の平面の2Dデカルト座標系である。サンプル110は、様々な視野角θを通じて回転され、放射線写真の系列が取得される。サンプル座標系rは、源座標系xに対して回転する。動的な要素μ_d(r,t)が単一の画像の取得中に変化しないことを仮定する。

Mxは検出器130の平面の2Dデカルト座標系である。ただし、Mは系の倍率である（R₂/(R₁+R₂)に等しい）。投影の近似が有効であることを仮定すると、検出器の座標Mx、視野角θ及び取得時間tで取得された画像強度Iは、以下により与えられる。

ただし、g(Mx,θ,t)は、以下により与えられる検出器130でのコントラストである。

ただし、Pは以下により定義される円錐ビームX線投影演算子を示す。

開示の再構成法は、画像I(Mx,θ,t)からの一連の再構成時点tにおいてN³ボクセルのグリッドのサンプル座標rで動的な（時間変化する）線減衰係数μ(r,t)を再構成する。

これを実現するために、静的な要素μ_s(r)及び動的な要素μ_d(r,t)が別々に再構成される。静的な要素μ_s(r)は、t<0の空き時間（すなわち、動的なプロセスがt=0で開始される前）に画像化されてもよい。この理由は、動的な要素μ_d(r,t)がゼロになるからである。これを行うために、πN/2の静的な投影画像I_s(Mx,θ)のセットが、全回転（complete revolution）に渡って4/Nラジアンで均等に分けられたπN/2の視野角θの単一のスキャン中に取得され、静的な要素μ_s(r)が、静的なコントラストg_s(x,θ)で標準の3DのCT再構成アルゴリズムを使用して再構成される。図１の円錐ビームの幾何学配置に適した１つの実装では、再構成アルゴリズムは、その計算効率のために選択されたFDK（Feldkamp-Davis-Kress）フィルターバック投影である。

ただし、F[]は、x₁に対する1Dフーリエ変換を示し、ξ₁はx1に対して双数（dual）であり、F^-1{}はξ₁に対する逆1Dフーリエ変換を示し、Bは以下のように画像セットh(Mx,θ,t)で定義されたFDKバック投影演算子を示す。

一般性を失うことなく、時間軸は逆になってもよい点に留意すべきである。時間t=0で終了する動的なプロセスを考えてもよく、プロセスが完了した後に減衰係数の静的な要素を再構成してもよい。更なる別法として、動的なプロセスが開始時間t=0と終了時間t=t_endとの間で生じる場合、初期及び最終の静的な画像の双方が再構成されてもよい。

他の画像化幾何学配置が使用される場合、式(6)の代わりに、これらの幾何学配置に適した他の既知の再構成法が使用されてもよい。

図２は、本発明の実施例に従って、図１の画像化幾何学配置を使用して、動的なプロセスがt=0で初期化された後の一連の離散取得時点t及び対応する離散視野角θで取得された投影画像I(Mx,θ,t)のセットから動的に変化するサンプルの線減衰係数の動的な要素μ_d(r,t)を再構成する方法200を示すフローチャートである。一般的に、離散視野角θは、サンプルの複数の全回転を表す。

方法200は、前述の仮定1及び3に整合した全ての解の空間に渡って、(i)解と測定データとの間のずれ、(ii)（仮定2のもとでの）解の時間微分の空間的局所化という品質指標を同時最適化することを試みる。この同時最適化は、従来のCS-CT再構成法の単純な拡張を通じて実現できない。CS-CT法は、典型的には解空間が凸状（convex）であることを仮定するが、前述の仮定1及び2に整合した解の空間は凸状ではない。

方法200の前処理ステップとして、単に動的な要素μ_d(r,t)による“動的な”コントラストg_d(Mx,θ,t)は、以下により画像I(Mx,θ,t)から取得される。

・式(4)の投影演算子Pを再構成された静的な要素μ_s(r)に適用し、“静的な”コントラストg_s(Mx,θ)を取得する。或いは、離散視野角θが、画像I_s(Mx,θ)が取得された視野角のサブセットである場合、視野角θでの測定された静的なコントラストg_s(Mx,θ)を直接利用する。

・取得されたコントラストg(Mx,θ,t)から（投影又は測定された）静的なコントラストg_s(Mx,θ)を減算し、動的なコントラストg_d(Mx,θ,t)を残す。

方法200は、動的なコントラストg_d(Mx,θ,t)で実行される。

動的な要素μ_d(r,t)での変化の時間スケールが、１つの全回転に渡ってπN/2の視野角で画像の単一のスキャンを取得する時間よりかなり大きい場合、各スキャンの中間における再構成時点tでの動的な要素μ_d(r,t)は、そのスキャン中に取得された動的なコントラストg_d(Mx,θ,t)から、例えば、フィルターバック投影（式(6)）により再構成されてもよい。これは、“総当たり（brute force）”手法と呼ばれることがあり、事実上、一連の独立した2D又は3D再構成を表す。しかし、関心のある用途では、動的な要素μ_d(r,t)の変化の時間スケールは、１つの全回転に渡ってπN/2の視野角で画像のスキャンを取得するのに要する時間と同程度又はそれ未満であるため、総当たり手法は適用できない。

その代わりに、方法200は、全回転毎にπN/2の画像I(Mx,θ,t)よりかなり少ないものから動的な要素μ_d(r,t)の安定した再構成を可能にするために、サンプルの線減衰係数μ(r,t)についての事前の認識を組み込む。

動的な要素μ_d(r,t)が再構成される一連の再構成時点は、有限の時間分解能Δにより分離される。各再構成時点Tにおいて、動的な要素μ_d(r,T)は、再構成時点Tの周囲の取得時点tに取得された連続する動的なコントラスト画像g_d(Mx,θ,t)の系列から再構成される。

１つの実装では、連接する再構成時点Tを分離する時間分解能Δは、サンプルの全回転に必要な時間の半分になるように選択される。この実装では、動的な要素μ_d(r,T)は、再構成時点Tを対称的に取り囲む取得時点tにおいて、サンプルの全回転に渡って取得された動的なコントラスト画像g_d(Mx,θ,t)の単一のスキャンから再構成される。すなわち、μ_d(r,T)の再構成に寄与するコントラスト画像g_d(Mx,θ,t)の半分が、Tの前の時点tに取得され、半分がTの後の時点tに取得される。従って、各コントラスト画像g_d(Mx,θ,t)は、動的な要素の２つの再構成に寄与する。１つはtの前の再構成時点Tであり、１つはtの後の再構成時点T+Δである。従って、この実装では、２つの再構成μ_d(r,T)及びμ_d(r,T+Δ)に寄与する連続するコントラスト画像g_d(Mx,θ,t)の連続するセットの間に50%の重複が存在する。この重複の存在のため、いずれか２つの連続する再構成された画像（例えば、μ_d(r,T)及びμ_d(r,T+Δ)）は、完全には独立していない。連続する再構成された画像のこの相互依存は、前述の“総当たり”手法に存在しない。連続する再構成された画像の相互依存は、再構成時点の間の変化が小さいという前述の過程2を実行するのに役立つ。他の実装において時間分解能Δを減少又は増加させることは、投影画像の連続する寄与するセットの間の重複の量を増加又は減少させる。

方法200は、CS及びDT法の双方の要素を組み込む。１つの実装では、方法200は、同時反復再構成法（SIRT：simultaneous iterative reconstruction technique）に基づく。SIRTは、限られたデータセットでのその既知の良好な性能のため選択される。方法200の他の実装では、アンダーサンプリング（under-sampling）の条件で適切にうまく動作する如何なる再構成及び再投影アルゴリズムが使用されてもよい。例は、反復フィルターバック投影（IFBP：iterated filtered backprojection）、SART及びARTを含むが、これらに限定されない。

図３Ａ及び３Ｂは、方法200が実施され得る汎用コンピュータシステム300の概略ブロック図を併せて形成している。

図３Ａに示すように、コンピュータシステム300は、コンピュータモジュール301と、キーボード302、マウスポインタデバイス303、スキャナ326、カメラ327及びマイクロフォン380のような入力デバイスと、プリンタ315、ディスプレイデバイス314及びスピーカ317を含む出力デバイスとにより形成される。外部変復調（Modem）送受信デバイス316は、接続321を介して通信ネットワーク320へ及び通信ネットワーク320から通信するためにコンピュータモジュール301により使用されてもよい。ネットワーク320は、インターネット又はプライベートWANのような広域ネットワーク（WAN）でもよい。接続321は、電話回線であり、モデム316は、従来の“ダイヤルアップ”モデムでもよい。或いは、接続321が高容量（例えば、ケーブル）接続である場合、モデム316はブロードバンドモデムでもよい。無線モデムがネットワーク320への無線接続のために使用されてもよい。

典型的には、コンピュータモジュール301は、少なくとも１つのプロセッサユニット305と、例えば、半導体ランダムアクセスメモリ（RAM）及び半導体読み取り専用メモリ（ROM）から形成されたメモリユニット306とを含む。モジュール301はまた、ビデオディスプレイ314、スピーカ317及びマイクロフォン380に結合するオーディオ・ビデオインタフェース307と、キーボード302及びマウス303のためのI/Oインタフェース313と、外部モデム316及びプリンタ315のためのインタフェース308とを含み、複数の入力／出力（I/O）インタフェースを含む。或る実装では、モデム316は、コンピュータモジュール301（例えば、インタフェース308）内に組み込まれてもよい。コンピュータモジュール301はまた、接続323を介してローカルエリアネットワーク（LAN）として知られるローカルコンピュータネットワーク322へのコンピュータシステム300の結合を可能にするローカルネットワークインタフェース311を有する。また、図示のように、ローカルネットワーク322はまた、典型的には“ファイヤウォール”デバイス又は同様の機能を有するデバイスを含む接続324を介して広域ネットワーク320に結合してもよい。インタフェース311は、Ethernet(TM)回路カード、Bluetooth(TM)無線構成又はIEEE803.11無線構成により形成されてもよい。

インタフェース308及び313は、直列及び並列接続の一方又は双方を提供してもよい。前者は、典型的にはUSB（Universal Serial Bus）標準に従って実装され、対応するUSBコネクタ（図示せず）を有する。記憶デバイス309が提供され、典型的には、ハードディスクドライブ（HDD：hard disk drive）310を含む。フロッピーディスクドライブ及び磁気テープドライブ（図示せず）のような他の記憶デバイスが使用されてもよい。外部の不揮発性データ源とインタフェース接続するために、典型的にリーダ312が提供される。例えば、光ディスク（例えば、CD-ROM、DVD）、USB-RAM及びフロッピーディスクのようなポータブルコンピュータの可読記憶デバイス325は、システム300への適切なデータ源として使用されてもよい。

典型的には、コンピュータモジュール301の構成要素305〜313は、相互接続バス304を介して、関連分野の当業者に既知のコンピュータシステム300の通常の動作モードを生じるように通信する。記載の構成が実施され得るコンピュータの例は、IBM-PC及びその互換機、Sun Sparcstation、Apple Mac(TM)又はこれらから進展したコンピュータシステムを含む。

以下に説明する方法は、コンピュータシステム300内で実行可能な１つ以上のソフトウェアアプリケーションプログラム333として、コンピュータシステム300を使用して実装されてもよい。特に、図３Ｂを参照すると、記載の方法のステップは、コンピュータシステム300内で実行されるソフトウェア333の命令331により実施される。ソフトウェア命令331は、１つ以上の特定のタスクをそれぞれ実行する１つ以上のコードモジュールとして形成されてもよい。ソフトウェアはまた、２つの別々の部分に分割されてもよく、第１の部分及び対応するコードモジュールは記載の方法を実行し、第２の部分及び対応するコードモジュールは第１の部分とユーザとの間のユーザインタフェースを管理する。

ソフトウェア333は、一般的にはコンピュータ可読媒体からコンピュータシステム300にロードされ、典型的には図３Ａに示すようなHDD310又はメモリ306に格納される。その後、ソフトウェア333は、コンピュータシステム300により実行されてもよい。或る場合には、アプリケーションプログラム333は、１つ以上の記憶媒体325にコード化されてユーザに供給されてもよく、メモリ310又は306に記憶する前に対応するリーダ312を介して読み取られてもよい。コンピュータ可読記憶媒体は、実行及び／又は処理のためにコンピュータシステム300に命令及び／又はデータを提供する際に関与するいずれかの非過渡的な有形の記憶媒体を示す。このような記憶媒体の例は、このようなデバイスがコンピュータモジュール301の内部にあっても外部にあっても、フロッピーディスク、磁気テープ、CD-ROM、DVD、ハードディスクドライブ、ROM又は集積回路、USBメモリ、光磁気ディスク、半導体メモリ、又はPCMCIAカード等のようなコンピュータ可読カードを含む。このようなソフトウェア又はコンピュータプログラムを記録したコンピュータ可読記憶媒体は、コンピュータプログラムプロダクトである。コンピュータモジュール301でのこのようなコンピュータプログラムプロダクトの使用は、動的な要素を再構成する装置を生じる。

或いは、ソフトウェア333は、ネットワーク320又は322からコンピュータシステム300により読み取られてもよく、他のコンピュータ可読媒体からコンピュータシステム300にロードされてもよい。コンピュータモジュール301へのソフトウェア、アプリケーションプログラム、命令及び／又はデータの提供にも関与し得る過渡的又は無形のコンピュータ可読伝送媒体の例は、無線又は赤外線伝送チャネルと、他のコンピュータ又はネットワーク接続されたデバイスへのネットワーク接続と、電子メール送信及びウェブサイトに記録された情報等を含むインターネット又はイントラネットとを含む。

アプリケーションプログラム333の第２の部分及び前述の対応するコードモジュールは、ディスプレイ314上に処理又は表示される１つ以上のグラフィカルユーザインタフェース（GUI：graphical user interface）を実装するために実行されてもよい。典型的にはキーボード302及びマウス303の操作を通じて、コンピュータシステム300及びアプリケーションのユーザは、GUIに関連するアプリケーションに制御コマンド及び／又は入力を提供するために、機能的に適応可能な方法でインタフェースを操作してもよい。スピーカ317を介した音声プロンプト出力を利用したオーディオインタフェース及びマイクロフォン380を介したユーザの音声のコマンド入力のように、機能的に適応可能なユーザインタフェースの他の形式も実装されてもよい。

図３Ｂは、プロセッサ305及び“メモリ”334の詳細な概略ブロック図である。メモリ334は、図３Ａのコンピュータモジュール301によりアクセス可能な全てのメモリデバイス（HDD310及び半導体メモリ306を含む）の論理集合を表す。

コンピュータモジュール301が最初に電源投入された場合、パワーオン・セルフテスト（POST：power-on self-test）プログラム350が実行する。POSTプログラム350は、典型的には半導体メモリ306のROM349に格納される。ROM349のようなハードウェアデバイスに永続的に格納されたプログラムは、場合によってファームウェアと呼ばれる。POSTプログラム350は、適切な機能を確保するために、コンピュータモジュール301内のハードウェアを検査し、典型的には、正確な動作のために、プロセッサ305と、メモリ（309、306）と、典型的にはROM349に格納されるBIOS（basic input-output systems software）モジュール351とを検査する。POSTプログラム350がうまく動作すると、BIOS351は、ハードディスクドライブ310を起動する。ハードディスクドライブ310の起動は、ハードディスクドライブ310に存在するブートストラップローダ352を、プロセッサ305を介して実行させる。これは、オペレーティングシステム353をRAMメモリ306にロードし、これにより、オペレーティングシステム353が動作を開始する。オペレーティングシステム353は、プロセッサ管理、メモリ管理、デバイス管理、記憶管理、ソフトウェアアプリケーションインタフェース及び一般ユーザインタフェースを含む様々なハイレベルの機能を実現するためにプロセッサ305により実行可能なシステムレベルのアプリケーションである。

オペレーティングシステム353は、コンピュータモジュール301で実行する各プロセス又はアプリケーションが他のプロセスに割り当てられたメモリと衝突することなく実行するのに十分なメモリを有することを確保するために、メモリ（309、306）を管理する。更に、システム300内で利用可能な異なる種類のメモリは、各プロセスが効率的に実行できるように適切に使用されなければならない。従って、集約されたメモリ334は、（特に言及しない限り）どのようにメモリの特定のセグメントが割り当てられるかを示すことを意図しているのではなく、コンピュータシステム300によりアクセス可能なメモリの一般的な図及びどのようにこのようなものが使用されるかを提供することを意図する。

プロセッサ305は、制御ユニット339と、論理演算装置（ALU：arithmetic logic unit）340と、場合によってはキャッシュメモリと呼ばれるローカル又は内部メモリ348とを含み、複数の機能モジュールを含む。典型的には、キャッシュメモリ348は、レジスタ部分に複数の記憶レジスタ344〜346を含む。１つ以上の内部バス341は、これらの機能モジュールを機能的に相互接続する。典型的には、プロセッサ305はまた、接続318を使用してシステムバス304を介して外部デバイスと通信するための１つ以上のインタフェース342を有する。

アプリケーションプログラム333は、条件分岐及びループ命令を含んでもよい命令331のシーケンスを含む。プログラム333はまた、プログラム333の実行の際に使用されるデータ332を含んでもよい。命令331及びデータ332は、それぞれメモリ位置328-330及び335-337に格納される。命令331とメモリ位置328-330との相対サイズに応じて、特定の命令は、メモリ位置330に示す命令により示されるように、単一のメモリに格納されてもよい。或いは、命令は、複数の部分に分割されてもよく、メモリ位置328-329に示す命令部分により示されるように、それぞれが別々のメモリ位置に格納される。

一般的に、プロセッサ305は、実行される命令セットを与えられる。プロセッサ305は、次の入力を待機し、次の入力に対して他の命令セットを実行することにより反応する。各入力は、入力デバイス302、303のうち１つ以上により生成されたデータ、１つ以上のネットワーク320、322を通じて外部源から受信したデータ、記憶デバイス306、309のうち１つから取得されたデータ、又は対応するリーダ312に挿入された記憶媒体325から取得されたデータを含み、複数の源のうち１つ以上から提供されてもよい。命令セットの実行は、或る場合には、データの出力を生じてもよい。実行はまた、メモリ334にデータ又は変数を格納すること含んでもよい。

開示の方法は、対応するメモリ位置355-358にてメモリ334に格納された入力変数354を使用する。開示の方法は、対応するメモリ位置362-365にてメモリ334に格納された出力変数361を生成する。中間変数は、メモリ位置359、360、366及び367に格納されてもよい。

プロセッサ305のレジスタ部分344-346、論理演算装置（ALU）340及び制御ユニット339は、プログラム333を構成する命令セットの命令毎に“取得、デコード及び実行”のサイクルを実行するために必要なマイクロ演算のシーケンスを実行するために共に動作する。それぞれの取得、デコード及び実行のサイクルは以下のものを有する。

(a)メモリ位置328から命令331を取得する又は読み取る取得動作
(b)どの命令が取得されたかを制御ユニット339が判定するデコード動作
(c)制御ユニット339及び／又はALU340が命令を実行する実行動作
その後、次の命令の更なる取得、デコード及び実行サイクルが実行されてもよい。同様に、制御ユニット339が値をメモリ位置332に格納する又は書き込む格納サイクルが実行されてもよい。

図２の方法の各ステップ又はサブプロセスは、プログラム333の１つ以上のセグメントに関連付けられ、プログラム333の記された部分について命令セットの命令毎に取得、デコード及び実行サイクルを実行するために共に動作するプロセッサ305のレジスタ部分334-337、ALU340及び制御ユニット339により実行される。

或いは、方法200は、効率的且つ並列に複数の小さい数学演算を実行する機能を用いて、グラフィックス処理ユニット（GPU：Graphics Processing Unit）に基づくコンピュータプラットフォーム又はコンピュータシステム300に類似した他のマルチプロセッサのコンピュータプラットフォームで実行されてもよい。

方法200はステップ210で始まり、動的な要素μ_d(r,t)の初期推定μ_d ⁽⁰⁾(r,t)が形成される。１つの実装では、動的な要素μ_d(r,t)の初期推定μ_d ⁽⁰⁾(r,t)は、ゼロに等しい。反復カウンタnもゼロに初期化される。

方法200はステップ220に進み、動的な要素μ_d(r,t)の現在の推定μ_d ⁽ⁿ⁾(r,t)が再構成時点に更新される。更新は、動的なコントラスト画像g_d(Mx,θ,t)を用いて更新された推定μ_d ⁽ⁿ⁾(r,t)の整合性を促すために、SIRT（又は他の適切なアルゴリズム、上記参照）の単一の反復を使用して行われる。更新を実行するために、以下に定義する動的なSIRT演算子Sが現在の推定μ_d ⁽ⁿ⁾(r,t)に適用される。

ただし、N(x,θ,t)は、再構成領域を通じた投影経路長であり、P及びBは、それぞれ式(4)及び(7)により定義される投影及び逆投影演算子である。式(9)の動的なSIRT演算子Sは、現在の推定μ_d ⁽ⁿ⁾(r,t)を投影し、動的なコントラストg_d(Mx,θ,t)からこれらの投影を減算し、投影経路長により差分を正規化し、正規化された差分を逆投影することにより、現在の推定μ_d ⁽ⁿ⁾(r,t)を更新する。

式(9)で使用するため、投影及び逆投影演算子P及びBの双方は、μ_d ⁽ⁿ⁾(r,t)の有限な時間分解能Δを考慮するため、時間軸に沿って補間を必要とする。投影ステップにおける補間は以下のように実行される。連続する動的な要素の推定μ_d ⁽ⁿ⁾(r,T)及びμ_d ⁽ⁿ⁾(r,T+Δ)からTとT+Δとの間の時点tにおいて投影された動的なコントラストg_d(Mx,θ,t)を計算するために、線形補間が使用される。

逆投影ステップで補間に使用される重みは、式(10)で与えられる投影重みの逆数である。

反復カウンタnもステップ220でインクリメントされる。

次のステップ230において、更新された推定μ_d ⁽ⁿ⁾(r,t)は、サンプルの動的な動作についての事前情報の仮定2に基づいて“変更訂正（change-corrected）”される。１つの実装では、動的な要素μ(r,t)の瞬間的な変化が空間的に局所化されるという仮定2の定式化に適して、演算子

が更新された推定μ_d ⁽ⁿ⁾(r,t)に適用される。ただし、

はソフト閾値演算子であり、

は時間に関する偏導関数である。演算子

は、１つの取得時点と次の時点との間の動的な要素の変化の空間的局所化を促す（CS用途における解の“希釈化（sparsification）”と同様である）。ソフト閾値演算子

の閾値は、データセットの予想信号対雑音比に比例するように選択される。

次のステップ235において、更新された推定μ_d ⁽ⁿ⁾(r,t)は、サンプルの静的及び動的な要素の相互作用についての事前情報の仮定3に基づいて、式(6)を使用して以前に再構成された静的な要素μ_s(r)を使用して“相互作用訂正（interaction-corrected）”される。ステップ235の１つの実装によれば、不浸透性の微孔性の静的な骨組を通じて非反応性流体の流れを画像化するときに適して、更新された推定μ_d ⁽ⁿ⁾(r,t)の空間支持領域が静的な要素μ_s(r)の空間支持領域の要素になることが仮定される。従って、更新された推定μ_d ⁽ⁿ⁾(r,t)は、その仮定の空間支持領域の外部でゼロに設定される。

方法200はステップ240に進み、方法が“よどみ（stagnation）”に到達したか否かを検査する。よどみは、現在の反復nのステップ220及び230での訂正後の更新された推定μ_d ⁽ⁿ⁾(r,t)の全体の変化が、典型的には、放射線写真の信号対雑音比に基づいて選択された或る値ε未満になるときに生じる。よどみが生じていない場合、方法200は他の反復のためにステップ220に戻る。よどみが生じた場合、方法200はステップ250に進む。

ステップ250において、訂正された推定μ_d ⁽ⁿ⁾(r,t)は、サンプルの動的な要素の物質組成についての事前情報の仮定1に基づいて“組成訂正（composition-corrected）”される。１つの実装では、動的な要素がバイナリ値になる（すなわち、２つの濃度値のみを使用して正確に表され得る）二相の非反応性非圧縮流体の流れの場合に適して、ステップ250は、訂正された推定μ_d ⁽ⁿ⁾(r,t)を“バイナリ化”する。ステップ250のこの実装は、以下のように定義されたバイナリ区分演算子Zを使用する。

バイナリ区分演算子Zは、訂正された推定μ_d ⁽ⁿ⁾(r,t)の絶対値が雑音閾値ε’未満になる場合に常に、訂正された推定μ_d ⁽ⁿ⁾(r,t)をゼロに設定する。その他の場合（すなわち、“ゼロでない領域”Ωでは、訂正された推定μ_d ⁽ⁿ⁾(r,t)は、ゼロでない領域Ωを通じて訂正された推定μ_d ⁽ⁿ⁾(r,t)の平均値を保持する値に設定される。バイナリ区分演算子Zは、ゼロでない領域Ωに渡って訂正された推定の値の事前の認識を必要としない。

３（又はn）相流体の流れの場合、ステップ250において、更に複雑な（すなわち、nレベルの）閾値演算（DT画像化から導かれる）が使用される。

バイナリ画像の空間は凸状でないため、各反復において組成訂正（ステップ250）を実行することは、方法200を誤った解に迅速に導く。

ステップ250の後に、ステップ260において、方法は、最近の更新、変更訂正、組成訂正ステップ220、230及び250の結合した効果が現在の推定μ_d ⁽ⁿ⁾(r,t)に有意な効果を有していないか否かを判定する。すなわち、ステップ260は、以下を検査する。

ただし、εはステップ240で使用された閾値である。正しい場合、方法200は終了する（ステップ270）。そうでない場合、方法200は他の反復のためにステップ220に戻る。

Australian National UniversityのX線マイクロCT機で実験データが収集された。“ビーズパック（beadpack）”サンプルは、ほぼ球状のAlSiO₂のビーズで詰め込まれた直径約1センチメートルのガラス管から形成されている。その結果の細孔の空間は、静的なサンプルを形成するために、コントラストについて0.5モルのヨウ化カリウムを投与されて、水で浸された。サンプルは、約68keVの特性ピークエネルギーで、2mmのSiO₂を通じてフィルタリングされたタングステンのターゲットから分岐する部分的に干渉性のX線で照射された。送信されたX線の強度は、画像センサとしてRoperのPI-SCX100:2048 X線カメラを使用して記録された。

720個の512×512ピクセルの放射線写真の全体の“静的な”スキャンは、１つの全回転に渡って均等に分けられた720の視野角で取得された。放射線写真毎の露光時間は1秒であった。静的なスキャンの完了時に、抽出ポンプがオンになり（時間t=0を定義する）、KIが投与された水が細孔の空間から流出した。第２の“動的な”放射線写真のセットは、水が流出されたときに収集された。動的なセットは、静的なスキャンと同様にそれぞれ1秒の露光時間で、均等間隔の視野角で、全回転毎に72個の放射線写真を有する。明らかに、静的なスキャンと比べて、“動的な”スキャンは、１つの完全なオーダーの大きさで角度についてアンダーサンプルされている。動的なデータ取得は、流体が完全に流出するまで続けられた。これは、30の全回転（約43分）を要した。

この例では、静的な線減衰係数μ_s(r)は、飽和したサンプルに対応し、動的な要素μ_d(r,t)は、流体が流出するときに形成する空間に対応する。前述の方法200を使用して再構成が実行された。動的なデータセットから再構成された3D画像の２つの2D視覚化が、図４に示されている。左の列は、T=15分21秒での再構成された3D画像の２つの表示を示しており、右の列は、左の列より１つの再構成時点だけ遅いT=15分59秒での再構成された3D画像の２つの表示を示している。各列の上の行は、異なる色でサンプルの静的な要素及び動的な要素を示しているが、下の行は動的な要素のみを示している。

図４の3D画像を比較すると、連続する再構成時点の間のビーズパックを通じたKIが投与された水の流出は、明確に見える。再構成された3D画像毎に38秒の実現された再構成時間分解能は、“総当たり”再構成法を使用して実現され得るもの（3D画像毎に15分）よりかなり優れている。

前述の構成は、特に、石油、地熱及び地滞留（geosequestration）業界に適用可能である。

前述のことは、本発明のいくつかの実施例のみを記載しており、本発明の範囲及び要旨を逸脱することなく、変更及び／又は変形が行われてもよい。実施例は例示的であり、限定的ではない。

Claims (20)

1. 一連の再構成時点で静的な要素と動的な要素とを有する動的なサンプルを表す多次元データセットを再構成する方法であって、
   前記動的なサンプルの複数の投影画像を取得するステップと、
   前記取得された投影画像から前記多次元データセットの前記静的な要素を再構成するステップと、
   前記動的なサンプルの更なる複数の投影画像を取得するステップと、
   前記動的なサンプルについての事前情報を使用して、前記更なる複数の取得された投影画像から各再構成時点での前記多次元データセットの前記動的な要素を再構成するステップと
   を有し、
   前記多次元データセットは、各再構成時点での前記静的な要素及び前記動的な要素の和である方法。
2. 前記第２の再構成は、
   前記更なる複数の投影画像を使用して、前記動的な要素の現在の推定を更新するステップと、
   前記動的なサンプルについての事前情報を使用して、前記動的な要素の前記更新された推定を訂正するステップと
   を有する、請求項１に記載の方法。
3. よどみが生じたか否かを判定するステップと、
   前記判定に応じて、前記訂正後の更新された推定の前記更新及び訂正ステップを繰り返すステップと
   を更に有する、請求項２に記載の方法。
4. 前記判定に応じて、前記動的な要素の前記訂正された推定を組成訂正するステップと、
   前記更新、訂正及び組成訂正が前記動的な要素の推定にかなりの効果を有しているか否かを判定するステップと
   を更に有する、請求項３に記載の方法。
5. 前記第２の判定に基づいて、前記更新、訂正及びよどみの判定を繰り返すステップを更に有する、請求項４に記載の方法。
6. 前記組成訂正は、前記訂正された推定の絶対値が雑音閾値未満になる場合には常に、前記訂正された推定をゼロに設定するステップを有する、請求項４に記載の方法。
7. 前記組成訂正は、前記訂正された推定の絶対値が雑音閾値以上になる場合には常に、前記訂正された推定を、前記訂正された推定の平均値を保持する値に設定するステップを有する、請求項６に記載の方法。
8. 前記訂正は、
   前記更なる複数の投影画像の取得時点の間に前記更新された推定の変化の空間的局所化を促すステップを有する、請求項２ないし７のうちいずれか１項に記載の方法。
9. 前記促すことは、前記多次元データセットの前記再構成された静的な要素の空間支持領域の補間である空間支持領域の外側で前記更新された推定をゼロに設定するステップを有する、請求項８に記載の方法。
10. 前記更新は、
    複数の視野角で前記動的な要素の現在の推定を投影し、複数の投影を形成するステップと、
    前記投影及び前記更なる複数の投影画像から差分画像を形成するステップと、
    前記サンプルを通じて投影経路長により前記差分画像を正規化するステップと、
    前記正規化された差分画像を逆投影するステップと、
    前記逆投影を前記動的な要素の前記現在の推定に加算するステップと
    を有する、請求項２ないし９のうちいずれか１項に記載の方法。
11. 前記形成は、
    前記更なる複数の投影画像から前記動的なサンプルの前記静的な要素の静的なコントラスト画像を減算し、前記動的なサンプルの動的なコントラスト画像を形成するステップと、
    前記動的なコントラスト画像から前記投影を減算するステップと
    を更に有する、請求項１０に記載の方法。
12. 前記動的なサンプルは、地層から抽出された地質標本である、請求項１ないし１１のうちいずれか１項に記載の方法。
13. 前記再構成された多次元データセットを使用して、前記地層の１つ以上の物理特性を判定するステップを更に有する、請求項１２に記載の方法。
14. 前記地層の前記判定された物理特性を使用して、前記地層から石油を抽出するステップを更に有する、請求項１３に記載の方法。
15. 複数の取得時点及び対応する視野角で取得された動的なサンプルの投影画像のセットから一連の再構成時点での動的なサンプルを表す一連の多次元画像を再構成する方法であって、各再構成時点で、
    前記視野角における前記再構成時点での前記多次元画像の現在の推定を投影し、複数の投影を形成するステップと、
    前記投影及び前記投影画像の系列から差分画像を形成するステップであり、前記系列は、前記再構成時点の周辺の取得時点で取得された連続する投影画像を有するステップと、
    前記動的なサンプルを通じて投影経路長により前記差分画像を正規化するステップと、
    前記正規化された差分画像を逆投影するステップと、
    前記再構成時点での前記多次元画像の現在の推定に前記逆投影を加算するステップと
    を有する方法。
16. 連続する再構成時点は、前記サンプルの全回転を構成する視野角で前記投影画像を取得する時間の半分に分けられる、請求項１５に記載の方法。
17. 前記系列は、前記サンプルの全回転を構成する視野角で取得された投影画像を有する、請求項１５に記載の方法。
18. 前記系列は、前記再構成時点を対称的に取り囲む取得時点において取得された投影画像を有する、請求項１５に記載の方法。
19. コンピュータデバイスに対して、一連の再構成時点で静的な要素と動的な要素とを有する動的なサンプルを表す多次元データセットを再構成する方法を実行させるように構成されたコンピュータプログラムであって、
    前記動的なサンプルの複数の投影画像を取得するコードと、
    前記投影画像から前記多次元データセットの前記静的な要素を再構成するコードと、
    前記動的なサンプルの更なる複数の投影画像を取得するコードと、
    前記動的なサンプルについての事前情報を使用して、前記更なる複数の取得された投影画像から各再構成時点での前記多次元データセットの前記動的な要素を再構成するコードと
    を有し、
    前記多次元データセットは、各再構成時点での前記静的な要素及び前記動的な要素の和であるコンピュータプログラム。
20. コンピュータデバイスに対して、複数の取得時点及び対応する視野角で取得された動的なサンプルの投影画像のセットから一連の再構成時点での動的なサンプルを表す一連の多次元画像を再構成する方法を実行させるように構成されたコンピュータプログラムであって、
    各再構成時点で、前記視野角における前記再構成時点での前記多次元画像の現在の推定を投影し、複数の投影を形成するコードと、
    前記投影及び前記投影画像の系列から差分画像を形成するコードであり、前記系列は、前記再構成時点の周辺の取得時点で取得された連続する投影画像を有するコードと、
    前記動的なサンプルを通じて投影経路長により前記差分画像を正規化するコードと、
    前記正規化された差分画像を逆投影するコードと、
    前記再構成時点での前記多次元画像の現在の推定に前記逆投影を加算するコードと
    を有するコンピュータプログラム。

Images