JP2014004359A - 反復式再構成の方法及び装置 - Google Patents

反復式再構成の方法及び装置 Download PDF

Info

Publication number
JP2014004359A
JP2014004359A JP2013127082A JP2013127082A JP2014004359A JP 2014004359 A JP2014004359 A JP 2014004359A JP 2013127082 A JP2013127082 A JP 2013127082A JP 2013127082 A JP2013127082 A JP 2013127082A JP 2014004359 A JP2014004359 A JP 2014004359A
Authority
JP
Japan
Prior art keywords
image
gradient
scaling factor
algorithm
voxel
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.)
Granted
Application number
JP2013127082A
Other languages
English (en)
Other versions
JP6212294B2 (ja
JP2014004359A5 (ja
Inventor
Baru Debasish
デバシッシュ・パル
Kristiaan Bernard De Man Bruno
ブルーノ・クリスティアン・バーナード・デー・マン
Fang Choi Jiang
ジャン・ファン・チョー
Yu Zhou
ツォー・ユ
Thibault Jean-Baptiste
ジャン−バプティスト・ティボルト
Fessler Jeffrey
ジェフリー・フェスラー
Fu Lin
リン・フー
Srivastava Somesh
ソメッシュ・スリヴァスタヴァ
Kim Donghwan
ドンファン・キム
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.)
General Electric Co
University of Michigan
Original Assignee
General Electric Co
University of Michigan
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 General Electric Co, University of Michigan filed Critical General Electric Co
Publication of JP2014004359A publication Critical patent/JP2014004359A/ja
Publication of JP2014004359A5 publication Critical patent/JP2014004359A5/ja
Application granted granted Critical
Publication of JP6212294B2 publication Critical patent/JP6212294B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • 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/408Dual energy
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/424Iterative

Landscapes

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

Abstract

【課題】診断画質を保ちつつCT走査の線量を有意に減少させるモデル・ベース反復式再構成(MBIR)アルゴリズムの計算時間及び計算資源を縮小する。
【解決手段】対象の画像を反復式で再構成する方法が提供される。この方法は、画像に関連する測定データを入手するステップと、画像を再構成するために同時型アルゴリズムを用いるステップとを含んでいる。画像を再構成するために同時型アルゴリズムを用いるステップは、ボクセル依存型であるスケーリング・ファクタを決定するステップと、画像を再構成するためにボクセル依存型スケーリング・ファクタを目的関数の勾配に適用するステップとを含んでいる。
【選択図】図1

Description

本書に開示される主題は一般的には、イメージング・システムに関し、さらに具体的には、反復式手法を用いて画像を再構成する方法及び装置に関する。
従来、画像は、計算機式断層写真法(CT)データからフィルタ補正逆投影(FBP)又は畳み込み逆投影(CBP)のような直接式再構成アルゴリズムを用いて再構成されている。近年、モデル・ベース反復式再構成(MBIR)アルゴリズムがCT画像を再構成するために出回るようになってきた。MBIRアルゴリズムの一つの利点は、これらのMBIRアルゴリズムがCTシステムから得られる測定をさらに正確にモデル化し得ることである。このことはマルチ・スライス検出器による螺旋CTシステムについて特にあてはまる。というのは、螺旋CTシステムは二次元(2D)再構成画像平面を斜めに通過する投影測定を発生し、取得されるデータは本質的に雑音を含むからである。これらの投影をさらに正確にモデル化することにより、MBIRアルゴリズムは、高品質(例えば高分解能)で、雑音が少なく、アーティファクトが少ない再構成を生成することができる。結果として、診断画質を保ちつつCT走査の線量を有意に減少させるツールとしてMBIRアルゴリズムを用いることができる。
しかしながら、MBIRの大きな問題は、再構成を完了するのに必要とされる計算時間及び計算資源である。MBIRアルゴリズムは典型的には、正確なシステム・モデル、統計学的雑音モデル、及び事前モデルを組み入れた目的関数を先ず形成することにより画像を再構成する。目的関数を利用し得るようにしておいて、次いで、目的関数を最小化する推定値を算出すること(典型的には反復式最適化アルゴリズムを用いて行なわれる)により画像を再構成する。かかる反復式最適化アルゴリズムの幾つかの実例として、反復式座標降下(ICD)法、期待値最大化(EM)法の諸変形、共役勾配(CG)法、及び順序付きサブセット(OS)法がある。しかしながら、MBIR目的関数及び関連する反復式解法は複雑であるため、幾つかの反復式アルゴリズムは、最終的な推定値に到達するために比較的多数回の反復を必要とし得る。結果として、MBIR目的関数を解く公知の反復式アルゴリズムは、画像を再構成するのに比較的長時間を必要とし得る。
一実施形態では、対象の画像を反復式で再構成する方法が提供される。この方法は、画像に関連する測定データを入手するステップと、画像を再構成するために同時型アルゴリズムを用いるステップとを含んでいる。画像を再構成するために同時型アルゴリズムを用いるステップは、ボクセル依存型であるスケーリング・ファクタを決定するステップと、画像を再構成するためにボクセル依存型スケーリング・ファクタを目的関数の勾配に適用するステップとを含んでいる。
もう一つの実施形態では、イメージング・システムが、検出器アレイと、該検出器アレイに結合されたコンピュータとを含んでいる。コンピュータは、対象の画像に関連する測定データを入手して、画像を再構成するために同時型アルゴリズムを用いるように構成されている。コンピュータは、ボクセル依存型であるスケーリング・ファクタを決定して、画像を再構成するためにボクセル依存型スケーリング・ファクタを目的関数の勾配に適用することにより、画像を再構成するために同時型アルゴリズムを用いるように構成されている。
もう一つの実施形態では、対象の画像を反復式で再構成する方法が提供される。この方法は、画像に関連する測定データを入手するステップと、画像を再構成するために同時型アルゴリズムを用いるステップとを含んでいる。同時型アルゴリズムは、正則化関数の勾配を算出することを含んでいる。画像を再構成するために同時型アルゴリズムを用いるステップは、複数のサンプル勾配値をルックアップ・テーブルに記憶させるステップと、ルックアップ・テーブルのサンプル勾配値を入手するステップと、サンプル勾配値を用いて近似的な実際の勾配値を算出するステップとを含んでいる。
もう一つの実施形態では、対象の画像を反復式で再構成する方法が提供される。この方法は、画像に関連する測定データを入手するステップと、画像を再構成するために同時型アルゴリズムを用いるステップとを含んでいる。同時型アルゴリズムは、正則化関数の勾配を算出することを含んでいる。画像を再構成するために同時型アルゴリズムを用いるステップは、画像の再構成の反復において実際の勾配を算出するステップと、実際の勾配に補償項を加えるステップと、画像の再構成の次回の反復のための勾配として実際の勾配及び補償項を用いるステップとを含んでいる。
様々な実施形態に従って対象の画像を再構成する方法の一例の流れ図である。 図1に示す方法を用いて生成される係数マップの一例の視覚的表現の図である。 図1の方法に関連する順序付きサブセット分離可能型放物面代理(OS−SPS)アルゴリズムの例の収束の例を示すグラフである。 図3のOS−SPSアルゴリズムの例を用いて形成される多様な画像の視覚的表現の図である。 様々な実施形態に従って対象の画像を再構成する方法の一例を示すもう一つの流れ図である。 螺旋幾何学的構成の一例及びこの螺旋幾何学的構成での再構成画像の一例の視覚的表現の図である。 図5の方法に関連するOS−SPSアルゴリズムの平均及び標準偏差の例を示すグラフである。 様々な実施形態に従って対象の画像を再構成するの方法の一例を示すもう一つの流れ図である。 様々な実施形態に従って対象の画像を再構成するの方法の一例を示すもう一つの流れ図である。 図9の方法に関連するOS−SPSアルゴリズムの一例を用いて形成される多様な画像の視覚的表現の図である。 図9の方法に関連するOS−SPSアルゴリズムの例に関連する収束速度を示すグラフである。 図9の方法に関連するOS−SPSアルゴリズムの例に関連する収束速度を示すもう一つのグラフである。 図9の方法に関連するOS−SPSアルゴリズムの例を用いて形成される多様な画像の視覚的表現の図である。 様々な実施形態に従って形成されるマルチ・モダリティ・イメージング・システムの一例の見取り図である。 図14に示すシステムのブロック概略図である。
以上の概要及び以下の様々な実施形態の詳細な説明は、添付図面と併せて読むとさらに十分に理解されよう。図面が様々な実施形態の機能ブロックの線図を示す範囲までにおいて、機能ブロックは必ずしもハードウェア回路の間の区分を示す訳ではない。従って、例えば、機能ブロックの1又は複数(例えばプロセッサ若しくはメモリ)が単体のハードウェア(例えば汎用信号プロセッサ若しくはランダム・アクセス・メモリのブロック、又はハードディスク等)として実装されてもよいし、多数のハードウェアとして実装されてもよい。同様に、プログラムは独立型プログラムであってもよいし、オペレーティング・システムのサブルーチンとして組み込まれていてもよいし、インストールされているソフトウェア・パッケージの機能等であってもよい。尚、様々な実施形態は図面に示されている構成及び手段に限定されないことを理解されたい。
本書には、同時型アルゴリズムを用いて画像を反復式で再構成する様々な実施形態が記載される。画像は複数の画像要素で構成される。本書に開示される目的のために、「画像要素」との用語は、画像空間配列の内部での対象の要素を指すものとする。例えば、画像要素は画像ピクセル又は画素を含むことができ、画像ピクセル又は画素は二次元(2D)又は三次元(3D)再構成での単一のボクセル又は体積要素に対応し得る。画像は、目的関数を最適化することにより再構成される。本書で用いられる目的関数又は費用関数は一般的には、撮像データを取得するのに用いられるイメージング・システムのモデル、このイメージング・システムに関連する雑音のモデル、及び再構成されている画像の形式のモデルを含んでいる。以下では、「最適化する」及び「最小化する」の各用語、並びに「目的関数」及び「費用関数」の各用語は互換的である。
様々な実施形態は、「適応」型のボクセル依存型スケーリング・ファクタを算出する新規のアプローチを提供する。ボクセル依存型スケーリング・ファクタは、画像を再構成するために同時型アルゴリズムの勾配に適用される。様々な実施形態は、取得(例えばスキャナ/リコン)の幾何学的構成に依存して画像に対する更新を緩和(又はスケーリング)する新規のアプローチを提供する。画像に対する更新は、勾配の適当なスケーリングを行なうこと、又は画像更新を算出するために勾配に適用されるステップ幅を修正することの何れかにより緩和(又はスケーリング)され得る。ステップ幅は典型的には、勾配に適用されるスカラーである。最適なステップ幅を算出すると最適な収束を齎すことができる。様々な実施形態は、基礎となる単調アルゴリズムの単調性に影響を与えることなく、正則化関数の勾配を算出する時間量を短縮する新規案を提供する。様々な実施形態の少なくとも一つの技術的効果は、計算時間を固定したままで過渡現象の低減及び画質の改善が得られることである。様々な実施形態の少なくとも一つの技術的効果は、従来の同時型アルゴリズムよりも少ない回数の反復で本来の費用関数の解にさらに高速に収束することである。例えば、様々な実施形態は、当技術分野で公開されている従来のアルゴリズムに必要とされる50回又は100回の反復に比較して、5回から10回の反復で許容可能な画質の収束尺度の類似値に収束し得る。
計算機式断層写真法(CT)撮像のためのモデル・ベース反復再構成(MBIR)アルゴリズムの開発は大きな注目を浴びている。というのも、MBIRアルゴリズムは雑音を低減して分解能を高めることにより画質を高めることが実証されているからである。MBIRアルゴリズムは典型的には、データ統計モデル、システムの幾何学的構成を記述するシステム・モデル、及び再構成画像の所望の特性を記述する事前モデルを含む費用関数を先ず形成することにより作動する。費用関数は典型的には複雑過ぎて測定から画像の直接的な再構成を1回で求めることができないため、画像は費用関数を何段階かで最小化することにより推定される。故に、かかるモデル・ベース費用関数を解くために反復最適化アルゴリズムが必要とされ、かかるアルゴリズムをMBIRアルゴリズムと呼ぶ。
μ∈RNを画像とし、y∈RMを測定された投影データとする。最大事後確率(MAP)推定問題は、データ対数尤度項Lと事前分布項Rとの和として定式化され得る。対数尤度項は、データの統計及びシステムの幾何学的構成に基づいて定式化される。事前分布項は、画像推定が何らかの確率密度関数(pdf)に属することを保証するものであり、βは事前分布スケーリング又は事前分布強度であってこの制約条件はβと共に適用する必要がある。
Φ(μ)=−ΦL(y,μ)+βΦR(μ) (1)
データ統計のモデルに依存して、対数尤度項Lの幾つかの選択肢が存在し得る。Lの一つのかかる選択肢は次の二次形態である。
ΦL(y,μ)=‖y−Aμ‖2 W (3)
式中、Wは、測定データの分散の近似的な逆行列として算出されるM×Mの対角行列であり、A∈RM×Nはシステムの幾何学的構成のモデルである。一連の反復最適化アルゴリズムの一方に、反復座標降下(ICD)法がある。ICDは、一度に1個の画像要素(「ボクセル」と呼ぶ)を、他の全てのボクセルを固定させたまま更新して、費用関数の勾配が1回1回のボクセル更新の後に更新されるような一連の一次元(1D)問題を作成する。1D費用関数の勾配及びヘッセ行列の更新を算出することは比較的簡単であり、良好な初期条件集合が与えられればICDは収束速度が比較的速い。ICDの一変種である非均一ICD(NH−ICD)は、幾つかのボクセルを他のボクセルよりも頻繁に更新して、さらに速い収束を実現する。ボクセルは一つのマップに基づいて選択され、このマップは全ての更新の履歴情報を用いて算出される。非均一ICDは、大きく更新されるボクセルはPWLS解から最も離隔しており、故にさらに頻繁な更新を必要とするとの観察に基づいている。しかしながら、ICD(及びこの各変種)は、高度並列化が可能なコンピュータ・プラットフォームであっても比較的効率よく比較的大規模に並列化することは困難である。
一連の反復最適化アルゴリズムの他方に、同時型アルゴリズムがある。本書で用いられる「同時型アルゴリズム」は、2個以上のボクセルを同時に更新するアルゴリズムである。また、「最適化」との用語は、画像を再構成する反復アルゴリズムの文脈での「更新」と同義である。MBIRの同時型アルゴリズムは、各々の段階でN次元の問題を解くことができる。同時型アルゴリズムは、比較的高度に並列化したコンピュータ・プラットフォームでの並列化に馴染み易い。かかるアルゴリズムの1回の反復はICDよりも高速に計算され得るが、同時型最適化アルゴリズムは通常、ICDよりも著しく多い回数の反復を収束に要する。故に、問題は、比較的速く収束し、且つ比較的効率のよい並列化に馴染み易い最適化アルゴリズムを見出すことである。
ペナルティ付き尤度関数Φの最小化因子(minimizer)を求める殆どの反復アルゴリズムが、尤度関数Φの勾配∇Φを利用する。多くのかかるアルゴリズムは以下の形態に書くことができる。
μn+1=μn+αnD(μn)∇Φ(μn) (4)
式中、αn>0は緩和係数(又はステップ幅)であり、N×N行列D(μ)はボクセル依存型スケーリング行列である。スケーリング行列D(μ)を本書では「スケーリング関数」及び/又は「スケーリング・ファクタ」とも呼ぶ。Dは最適化アルゴリズムに依存して対角であってもなくてもよい。従来の勾配降下アルゴリズムでは、Dは単位行列であり、αnを適正に選択することにより収束が保証され、αnは直線探索と呼ばれる1D最適化を用いて算出される。しかしながら、勾配降下は比較的収束が遅いため、実用目的には不満がある。共役勾配(CG)は、よりよい探索方向を求める観点では勾配降下に勝る。CGアルゴリズムを前処理すると、反復回数の関数として収束速度の高速化が得られる。理想的な前処理行列は費用関数のヘッセ行列の逆行列であり、この場合にはこの方法はニュートン法となり、二次の問題については単一段階で解かれ得る。しかしながら、断層写真法問題では、ヘッセ行列は空間不変であり、比較的大きく、反転が比較的困難である。費用関数のヘッセ行列のために明示型の対角行列及び循環行列による近似が提案されている。また、暗黙型前処理行列を導くこともできる。
また、現在の反復μnを用いることによりμn+1を算出するために、最適化転換(optimization transfer)の原理を用いて、最大化が比較的容易な代理関数φ(・;μn)を構築して問題を解くことができる。代理関数は以下の特性を有するように構築される。
Φ(μ)≦φ(μ,μn
Φ(μn)=φ(μ,μn)|μ=μn
∇Φ(μ)|μ=μn=∇10φ(μ,μn)|μ=μn (5)
分離可能型放物面代理(SPS)アルゴリズム、交互最小化アルゴリズム、及び期待値最大化アルゴリズムのような最適化転換アルゴリズムは代理関数を用いる。
「ブロック反復」法(「順序付きサブセット」法とも呼ぶ)が提案されている。かかるアルゴリズムは、適正なブロック依存型スケーリングによって多くの同時型アルゴリズムの収束を著しく高速化することができる。順序付きサブセット・アルゴリズムは、前方投影演算のレンジ分解を用いて基礎となるアルゴリズムを高速化する。各々のサブセットはデータのサブセットから成っており、データのサブセットを用いて全勾配を近似することにより高速化が達成される。これにより、大雑把にサブセット数倍だけ高速化する。しかしながら、順序付きサブセットで高速化されるアルゴリズムは単調型でなく、望ましくないリミット・サイクルに入る。リミット・サイクルの大きさはサブセット数に依存し、サブセットが大きいほどリミット・サイクルが大きくなる。リミット・サイクルは大きいほど望ましくなく、アルゴリズムが終結したときに実際の望ましい解から著しく逸脱する場合がある。
同時法及び「ブロック反復」法は典型的には、高度並列化が可能であるが、NH−ICDにおいて提案されているように履歴情報に基づいて更新に影響を与えるという融通性を欠く場合がある。かかるアルゴリズムは通常、低周波収束が比較的速く、高周波収束が比較的遅い。NH−ICDは、収束から離隔したボクセルほど頻繁に更新することにより、ベースラインICDアルゴリズムの収束を著しく高速化することが銘記されている。図1に示す方法100では、更新のために「適応」型のボクセル依存型スケーリング・ファクタD(μ)(式4)を算出し、従来の勾配に基づく最適化アルゴリズムよりも速く収束する新規のアプローチについて議論する。公知の再構成方法では、項「再スケーリング」Dを用いて「ブロック反復」アルゴリズムの収束を高速化しているが、再スケーリング・ファクタは幾何学的構成によって決まる固定的なファクタであり、反復番号と共に変化する訳ではなく、データに適応しない。本書に記載し且つ/又は図示する方法では、提案されるスケーリング・ファクタはデータに適応し、反復回数と共に変化し、履歴更新、取得幾何学的構成、データ値、画像特徴、及び/又は同様のもののようなファクタに基づいて算出され得る。
螺旋CTのような幾何学的構成では、ボクセルは走査軌跡によって一様にサンプリングされる訳ではない。典型的には、利用者によって予め設定された三次元(3D)着目領域(ROI)が再構成される。ビューの数はROIによって決まり、ROIの外部のボクセルの方がROIの内部のボクセルよりも少ない投影ビューによってサンプリングされ得る。収束過程の間に、これらのボクセルは顕著な過渡的アーティファクトを呈し、中間推定の画質を損なう場合がある。図5に示す方法200では、取得幾何学的構成に依存して、更新μn+1−μnを緩和(又はスケーリング)する(式4を参照されたい)新規のアプローチについて記載し図示する。これらの更新は、勾配∇Φの適当なスケーリングを行なうこと又はステップ幅αnを修正することの何れかによって緩和(又はスケーリング)され得る。方法200は過渡現象を低減し、画質を高めることができる。文献には「ブロック反復」アルゴリズムの収束を保証する緩和係数αnについて様々な方式が提案されている。しかしながら、文献の文脈でのスケーリング・ファクタは反復回数に依存する一連のスカラーである。本書に記載されるスケーリング・ファクタは空間変化する値であって、取得の幾何学的構成、履歴更新、データ値、画像特徴、及び/又は同様のもののような要因に依存し得る。
式1及び式4を参照すると、MBIR費用関数のための反復最適化アルゴリズムは、各回の反復で正則化項(regularizer)R(μ)の勾配を算出することを必要とする。正則化関数は典型的には、羃、対数、及び三角関数から成る。従って、正則化項の勾配を算出するのに必要とされる時間はかなり長くなる場合がある。図8及び図9にそれぞれ示す方法300及び方法400では、基礎となる単調アルゴリズムの単調性又はアルゴリズムの最適化の収束速度に影響を与えずに、正則化項の勾配を算出する時間を短縮する新規案を提示する。
図1に示す方法100では、式4の対角行列D(μ)をスケーリングするのに用いられる「適応」型のボクセル依存型マップを効率よく構築するために用いられて、MBIR費用関数のための反復最適化アルゴリズムの高速化をさらに齎し得る様々なアプローチについて議論する。マップの適応性を高めると共にマップをさらに効率よく構築する幾つかの案について議論する。ボクセル依存型マップを用いて対角行列D(μ)を再スケーリングする一例を、順序付きサブセットSPS(OS−SPS)アルゴリズムを用いて示す。
図5に示す方法200では、更新(μn+1−μn)を適正に緩和(又はスケーリング)して、過渡現象の低減及び画質の向上を齎す手法について記載する。取得幾何学的構成を用いて適当な緩和(又はスケーリング)ファクタを算出する一例を、OS−SPSアルゴリズムについて記載する。
図1は、対象の画像を反復式で再構成する方法の一例100の流れ図である。画像に関連する測定データが例示的なイメージング・システムを用いて取得される。ブロック102では、方法100は、測定データを入手するステップを含んでいる。ブロック104では、画像を再構成するために同時型アルゴリズムを用いる。画像を再構成するために同時型アルゴリズムを用いるブロック104のステップは、ブロック106においてボクセル依存型であるスケーリング・ファクタD(μ)を決定することを含んでいる。方法100のこの実施形態の例では、スケーリング・ファクタは、ブロック108において係数マップを生成することにより、ブロック106において決定される。Thibault他によって、解から離隔したボクセルほど頻繁に更新することにより収束が高速化し得ることが特記されている。故に、かかるボクセルを識別する必要があり、また「係数マップ」はかかるボクセルから成るべきである。NH−ICDでは、ボクセルの更新順序を選択するように係数マップを直接組み入れることができる。同時型アルゴリズムでは、係数マップを組み入れるように対角スケーリング行列D(μ)を修正することができる。アルゴリズムの単調性が確実に変化しないように注意を払うべきである。
かかる係数マップの一つの自然な選択肢は、少なくとも2回の反復又は少なくとも2回の下位反復についての差である。少なくとも2回の反復又は下位反復は連続していてもいなくてもよい。一実施形態では、係数マップは2回の連続した反復の差であって、下式によって与えられる。
j n=|μj n+1−μj n1 (6)
式6の差は、最も大きく変化し故に解から離隔している確率の高いボクセルを強調するものとなる。この差は反復回数と共に変化する。故に、係数マップは反復回数の関数となり、データに適応し得る。係数マップuj nは、式6の差のヒストグラムに基づいて変換され得る。この変換は、離散型であっても連続型であってもよい。理想的なマップは局所的に滑らかであり、故に最終的なマップは中間マップから低域通過フィルタ処理によって得ることができる。しかしながら、係数マップは1回の完全な反復が終わるまで利用可能にならないので、初期係数マップを初期画像から生成することができる。初期画像は典型的には、フィルタ補正逆投影(FBP)アルゴリズムを用いて再構成されるが、選択肢はFBPに限らない。初期画像にエッジ検出演算若しくは強度検出演算又は両方を施して用いて初期マップを生成することができ、この初期マップはまた幾何学的構成にも依存し得る。初期画像はかなりの低周波数内容を有すると想定すると、エッジ検出演算及び強度検出演算は、高周波収束を改善する妥当な選択肢であると言える。
係数マップは、画像容積についてスケーリングを適応式で選択する堅実な枠組みを提供する。係数マップを構築する判断は、限定しないが初期画像に関する事前知識、取得幾何学的構成に関する情報、現在の画像の特徴のような画像特徴(限定しないがエッジ、均一領域、及び/若しくは同様のもの等)、並びに/又は同様のもののような幾つもの要因を用いて強化され得る。例えば、FBPアルゴリズムは走査視野の外部の走査対象の切り捨てのためアーティファクトを呈し、この情報を利用して、切り捨てられた領域での更新をさらにスケーリングすることができる。もう一つの例は、同時型アルゴリズムは高周波収束が比較的遅いとの事実によって促されるものである。故に、あらゆる反復回での係数マップunが、μnの高域通過フィルタ処理済み画像から得られるマップによって強化され得る。
係数マップuj 4:uj nの一例を図2に示しており、同図では、係数マップがマップのヒストグラムを用いて離散化されている。最大の5%のボクセルに16の値を割り当て、次の10%、20%、及び40%のボクセルにそれぞれ8、4、及び2の値を割り当てている。ボクセルの残りには1の値を割り当てた後に、低域通過フィルタ処理を施している。通例のSPS更新は全てのボクセルを同等に更新するが、適応スケーリングされたSPS更新は、図で明るいボクセルほど大きい更新を施している。
図1へ戻り、係数マップの生成の後に、係数マップを用いて対角スケーリング行列D(μ)に作用させる。従って、方法100はブロック110において、ボクセル依存型スケーリング・ファクタD(μ)を算出するために係数マップを用いるステップを含んでいる。一旦ボクセル依存型スケーリング・ファクタD(μ)が決定されたら、方法100はブロック112において、画像を再構成するためにボクセル依存型スケーリング・ファクタD(μ)を目的関数の勾配∇Φに適用するステップを含んでいる。
係数マップは、本来のアルゴリズムの単調性に影響を与えずにボクセル依存型スケーリング・ファクタD(μ)に作用するものであるべきである。本来のアルゴリズムが単調でないとき、例えば「ブロック反復」アルゴリズムでは気にしなくてよい。無視すべきでない実際的な観点は、具現化形態において対角スケーリング行列をスケーリングするために係数マップを組み入れる効率である。以下、幾つかの可能な案について記載する。
対角スケーリング行列D(μ)をスケーリングするために係数マップを利用するという案は、限定しないが暗黙型又は明示型前処理行列によるCGアルゴリズムのような同時型アルゴリズムの殆どに用いることができる。ここで、方法100の一実施形態を例示するOS−SPSアルゴリズムの一例について記載する。OS−SPSアルゴリズムは、ボクセル依存型係数マップを用いて対角行列D(μ)を再スケーリングする一例である。
ペナルティ付き重み付き最小自乗費用関数は、下式のように書かれる。
L(y,μ)=−(1/2)Σi=1 Mi(yi−[Aμ]i2
βΣr=1 KΨr([Cμ]r),K=13×N,C∈RK×N
=(1/2)Σi=1 Mi([Aμ]i)+βΣr=1 KΨr([Cμ]r) (7)
式中、yはM個の要素を有する投影データであり、wは統計学的重み付け行列であり、Aはシステム行列であり、μはN個の要素を有する画像であり、Cは有限差分行列であり、Ψは有向係数を含む任意のマルコフ確率場(MRF)事前確率である。分離可能型代理アルゴリズムは、De Pierroによって提案された凸分解補題によって促されたものであり、代理関数のヘッセ行列の比較的容易な反転を可能にすると共に、非負制約を組み入れることが比較的容易である。以下の変換を用いてデータ・フィット項及び正則化項の代理式を導くことができる。
[Aμ]i=Σj=1 Nijμj
=Σj=1 Nπij n[(aij/πij n)(μj−μj n)+[Aμni
[Cμ]r=Σj=1 Nrjμj
=Σj=1 Nκrj n[(crj/κrj n)(μj−μj n)+[Cμ]r] (8)
式中、Σj=1 Nπij n=1及びΣj=1 Nκrj n=1である。式8を用いると、代理式を下式のように構築することができる。
i([Aμ]i)≦Qi(μ,μn
=Σj=1 Nπij ni([(aij/πij n)(μj−μj n)+[Aμni])
Ψr([Cμ]r)≦Rr(μ,μn
=Σj=1 Nκrj nΨr([(crj/κrj n)(μj−μj n)+[Cμnr])
(9)
正則化項(regularizer)の代理式についてπij=aij/Σjij及び最大曲率を用いるとSPSアルゴリズムの一形態が導かれ、下式のように書くことができる。
μj n+1=μj n
(Σj=1 M∇qi([Aμni)+Σr=1 K∇Ψr([Cμnr))
/(dj Q,n+dj R,n
j Q,n=ΣiijiΣjij
j R,n=Σrrj2Ψr(0)Σjrj (10)
式10では、対角スケーリング行列は、代理費用関数(式9)のヘッセ行列の逆行列である。このヘッセ行列は対角行列であって、故に比較的容易に反転させることができる。このスケーリング行列を選択すると、スケーリング行列はデータ統計及び幾何学的構成の両方にのみ依存するものとなる。故に、スケーリング行列は固定され、反復番号と共に変化することはない。
πij及びκrjを修正することにより非一様マップujをスケーリング行列に組み入れることができる。πij=aijj/Σjijj及びκrj=crjj/Σjrjjを用いると、非一様更新によるSPS(NU−SPS)の更新を下式のように書くことができる。
μj n+1=μj n
(Σi=1 M∇qi([Aμni)+Σr=1 K∇Ψr([Cμnr))
/(dj Q,n+dj R,n
j Q,n=(ΣiijiΣjijj)/uj
j R,n=(Σrrj2Ψr(0)Σjrjj)/uj (11)
jを含めると、反復回数に依存し、データに適応し得る適応項dj・,nが導かれる。式11は、一つの完全な前方投影及び逆投影を用い、係数マップuj nを用いてdj・,nを算出している。図3(a)は、患者データについて収束済み画像からの自乗平均平方根(RMS)差(HU単位)を反復回数に対してプロットしており、対角行列の適応スケーリングによって収束速度が約3倍高速化することを示している。また図4は、新たなスケーリング行列による再構成画像が、新たなスケーリングを行なわない場合よりも収束済みの画像に近いことを確認している。図4(a)は初期FBP画像を示している。図4(b)は収束済みの画像を示している。図4(c)及び図4(d)は、(c)通常のOS−SPS(82個のサブセット)及び(d)適応スケーリングされたOS−SPS(82個のサブセット)による20回反復の後の再構成画像を示している。また、図3(b)は、高周波成分の収束速度が3倍を上回って高速化していることを示している。(高周波項の収束は、高コントラスト・ワイヤの半値幅(FWHM)を算出することにより測定されている。)
無視すべきでない実際的な観点は、対角スケーリング行列をスケーリングするために係数マップを組み入れる効率である。以上の案は具現化形態毎に特定的であり、アルゴリズム間で異なり得る。計算量を小さくする幾つかの可能な案は以下の通りである。
1.係数マップを算出する頻度を減少させ、各回の反復でなく数回の反復(例えばn回の反復)毎に算出することができる。
2.係数マップをダウンサンプリングして前方投影及び逆投影の時間を短縮することができる。係数マップは局所的に滑らかであり、故にダウンサンプリングは収束に大きな影響がないと言える。
3.式11の更新は勾配∇qiを算出することを必要とし、前方投影及び逆投影を必要とする。πijを算出するための前方投影及び逆投影を勾配計算と同時に行なうことができる。これにより、対角スケーリング行列D(μ)を算出するというオーバヘッドを減少させることができる。例えばOSアルゴリズムでは、勾配は部分的な前方投影及び逆投影を用いて算出される。この場合には、新たなスケーリング行列を勾配と共に算出することができるが、この完全な計算は、サブセットの全てではないとしても少なくとも幾つかを通じてのサイクル計算を必要とし得る。
図5は、対象の画像を反復式で再構成する方法のもう一つの例200の流れ図である。マルチ・スライス螺旋CTのような幾つかの取得幾何学的構成では、画像容積は非一様にサンプリングされる。図6(b)に示すように、三次元着目領域(3D−ROI)の内部のボクセルは、ROIの外部のボクセルに比較して多い投影ビューによってサンプリングされる。これらアンダーサンプリングされたボクセルのため、アルゴリズムの収束が著しく遅くなる場合があり、現在の画像に対する更新(式4)が適当にスケーリングされていなければ過渡的なアーティファクトを招き得る。
図5へ戻り、画像に関連する測定データが例示的なイメージング・システムを用いて取得される。ブロック202では、方法200は、測定データを入手するステップを含んでいる。測定データは、3D取得幾何学的構成において非一様にサンプリングされている。ブロック204では、画像を再構成するために同時型アルゴリズムを用いる。画像を再構成するために同時型アルゴリズムを用いるブロック204のステップは、ブロック206においてボクセル依存型であるスケーリング・ファクタD(μ)を決定することを含んでいる。方法200のこの実施形態の例では、ボクセル依存型スケーリング・ファクタは、取得幾何学的構成に基づいて構築される。従って、ブロック208では、方法200は取得幾何学的構成に基づいて直接ボクセル依存型スケーリング・ファクタを算出することを含んでいる。
一旦ボクセル依存型スケーリング・ファクタD(μ)が決定されたら、方法200はブロック210において、画像を再構成するためにボクセル依存型スケーリング・ファクタD(μ)を目的関数の勾配∇Φ(式4を参照されたい)に適用するステップを含んでいる。安定性のために勾配∇Φをスケーリングしてもよいしステップ幅αnを修正してもよい。例えば、幾つかの実施形態では、ボクセル依存型スケーリング・ファクタを勾配に適用するブロック210のステップは、ボクセル依存型スケーリング・ファクタを用いて勾配をスケーリングすることを含んでいる。このスケーリングは、ATWAを用いて得ることができる(式15を参照されたい)。スケーリングを一つ一つのボクセルに対するサブセットの寄与と関係付けて、更新の適当なスケーリングを導くことができる。他の実施形態では、ボクセル依存型スケーリング・ファクタを勾配に適用するブロック210のステップは、画像データの各々のスライスについて緩和係数を算出することを含んでいる。画像容積のサンプリングは非一様であり、ROIの内部のスライスはROIの外部のスライスに比較して多くのビューによってサンプリングされている。ステップ幅(又は緩和係数)は通常スカラーであるが、各々のスライス毎にベクトルとなるように修正することもできる。これによりROIの外部のスライスをROIの内部のスライスよりも緩和させることが可能になり得る。しかしながら、これによりαを求める探索空間の次元が増加し、計算時間が増大する場合がある。
アンダーサンプリングされたボクセルについて更新を緩和(又はスケーリング)するという案は、全ての同時型アルゴリズムに適用可能である。ここで、方法200の一実施形態を示すOS−SPSアルゴリズムの一例について記載する。OS−SPSアルゴリズムは、スケーリング・ファクタの構築、及び勾配∇Φをスケーリングするためのスケーリング・ファクタの利用の一例である。
式10を用いて、SPSアルゴリズムについての更新を下式のように書くことができる。
μj n+1=μj n
[Σiiji(Σjijμj n−yi)+βΣr∇Φr(Σjrjμj n)]
/dj (12)
式中、djは、代理費用関数のヘッセ行列から決定されるスケーリング関数である(式4を参照されたい)。djは、式10に記述されているデータ・フィット項及び正則化項から成る。さて、OS−SPSにおいて、勾配はビューのサブセットを用いて算出され、1回の完全な反復は全S個のサブセットを通じたサイクル計算に対応する。各回の反復は多数回の下位反復を含んでおり、各回の下位反復は、データのサブセットを用いて正則化項の勾配及びデータ・フィット項を算出することを必要とする。OS−SPS更新は、下式のように書くことができる。
μj n,s+1=μj n,s
[γj n,sΣiisjis(Σjisjμj n,s−yis)+
βΣr∇Φr(Σjrjμj n)]
/dj (13)
式中、sはサブセット番号を示し、nは反復番号を示す。適正にスケーリングされれば、ビューのサブセットを用いて算出される勾配は、全てのビューを用いて算出される勾配に一致し得る。これにより、本来のSPSアルゴリズムに比較して収束がS倍高速化される。しかしながら、この理想的なスケーリング・ファクタ
γj n,s=[Σisisjis(Σjisjμj n,s−yis)]
/[Σiiji(Σjijμj n,s−yi)] (14)
を算出することは実際的でない。OS−SPSの従来の実際的な具現化形態は、サブセット数に対応する定数スケーリング・ファクタγj n,s=Sを用いるものである。しかしながら、この選択肢は、理想的なスケーリング・ファクタγj n,sが著しく非一様であるときには良好に作用しない。スケーリングのこの非一様性のため、アルゴリズムが望ましくない顕著な過渡的アーティファクトを有するものになる場合がある。この非一様性の理由の一つは取得幾何学的構成であり得る。取得幾何学的構成が画像容積の非一様サンプリングを招く場合がある。勾配に対するボクセル依存型スケーリングは、取得幾何学的構成に基づいて展開され得る。非一様マップの一つのかかる選択肢は、
γj=Σs=1 SI{Σisisjis(Σjisj)>0} (15)
である。式中、I[B]はBが真であるときには1であり、Bが偽であるときには0である。式15では、jによって番号付けされた所与のボクセルがsによって番号付けされたサブセットに属するビューからの寄与を受けるときに、Bは真となる。非一様マップ及び分母dj=ΣiijiΣjij=Σs=1 SΣisisjis(Σjisj)は同時に計算され得る。
図6(b)は、適正にスケーリングされたOSアルゴリズムが、ROIの外部で相対的に急速には収束しない通常のOSアルゴリズムに対して高められた画質を提供することを示す。ROIの外部に現われた不安定性は、ROIの内部の画質に影響を与える(図6(b)の雑音標準偏差によって分かる)が、提案された緩和(又はスケーリング)によって減少している。図6(a)は、螺旋幾何学的構成の図を示している。斜線の領域は、ROIのボクセル及びROIの外部のボクセルの両方から寄与を受ける測定データを取得する検出器を示す。図6(b)は、螺旋幾何学的構成におけるOS−SPSアルゴリズムでの勾配スケーリングの効果を示す。各々の画像は、通常のスケーリング・アプローチ及び提案しているスケーリング・アプローチを用いて、328個のサブセットによるOSアルゴリズムを30回反復して実行した後に再構成されている。比較のために一様な領域(枠で囲まれている)の標準偏差σを算出している。(収束アルゴリズムの数回の反復を参照として示す。)
図7は、ROIの内部での通常のOSアルゴリズムはROIの外部での不安定性のため劣化するが、適正にスケーリングされたOSはROIの内部で堅牢であることをさらに証明している。図7は、ROIの末端スライスの一様な領域の内部での平均及び標準偏差を反復回数に対して示しており、提案しているOSアプローチ(328個のサブセット)に比較して通常のOSアプローチの不安定性を示す。
図8は、対象の画像を反復式で再構成する方法のもう一つの例300の流れ図である。図8は、正則化項勾配計算の実施形態の一例を示す。MBIR費用関数のための反復最適化アルゴリズムは、現在の画像推定値を更新するために正則化関数の勾配を算出することを必要とする(式4を参照されたい)。正則化項の勾配を算出する時間は、正則化項ポテンシャル関数Φに依存してかなりのものとなり得る。最も一般的には、このポテンシャル関数は羃、対数、及び三角関数を含んでおり、かなりの計算を必要とする。かかるポテンシャル関数の例として、限定しないが下記のものがある。
ΦR(μ)=log(1+|μ|2) (16)
ΦR(μ)=log(cosh(|μ|)) (17)
ΦR(μ)=[(1/2)|μ|p]/[1+|μ/δ|(p-q)] (18)
正則化関数の勾配を算出するのに必要とされる費用を減少させるために、正則化項勾配∇ΦR(μ)の近似を、方法300に記載するようにルックアップ・テーブルのサンプル勾配値に基づいて用いることができる。代替的には、図9の方法400に記載するように、正則化関数についての代理式を構築することができ、この場合には代理式の勾配の方が実際の正則化関数の勾配よりも効率よく算出され得る。
ここで方法300を参照すると、画像に関連する測定データが例示的なイメージング・システムを用いて取得される。ブロック302では、方法300は測定データを入手するステップを含んでいる。ブロック304では、画像を再構成するために同時型アルゴリズムを用いる。画像を再構成するために同時型アルゴリズムを用いるブロック304のステップは、ブロック306において複数のサンプル勾配値をルックアップ・テーブルに記憶させるステップを含んでいる。選択随意で、ブロック306の記憶させるステップは、固定されたサンプリング周波数で勾配値をサンプリングすることを含んでいる。ブロック306でサンプル勾配値をルックアップ・テーブルに記憶させるための一つのアプローチは、∇ΦR(μ)の値をサンプリングする、すなわちk=0,…,Kについてdk=∇ΦR(μk)をテーブルとして構成するものである。Δがサンプリング周波数を表わす正の整数である場合に、点μk=kΔであるならば、テーブルについて以下のインデクス生成関数を用いることができる。(式中、
は床関数を示す。)
ブロック308では、方法300は、ルックアップ・テーブルのサンプル勾配値を入手するステップを含んでいる。次いで、ブロック310では、ルックアップ・テーブルからのサンプル勾配値を用いて近似的な実際の勾配値を算出する。ブロック310での近似的な実際の勾配値は、限定しないが線形補間方式、二変数三次補間方式、及び/又は同様のもののような任意の補間方式を用いて算出され得る。一実施形態では、近似的な実際の勾配値を算出するブロック310のステップは、再構成の各回の反復毎に近似的な実際の勾配値を算出することを含んでいる。
図9は、対象の画像を反復式で再構成する方法のもう一つの例400の流れ図である。図9は正則化項勾配計算の実施形態の一例を示しており、上で概略的に述べたように正則化関数について代理式を構築する。画像に関連する測定データが例示的なイメージング・システムを用いて取得される。ブロック402では、方法400は測定データを入手するステップを含んでいる。ブロック404では、画像を再構成するために同時型アルゴリズムが用いられる。
画像を再構成するために同時型アルゴリズムを用いるブロック404のステップは、ブロック406において、画像の再構成の反復において実際の勾配を算出するステップを含んでいる。次いでブロック408では、実際の勾配に補償項を加えることにより正則化関数の代理式を構築する。補償項は比較的単純な数学的演算から成っていてよい。例えば、幾つかの実施形態では、補償項は加算、減算、除算、及び/又は乗算の少なくとも一つを含むのみである。
次いでブロック410では、実際の勾配及び補償項を、画像の再構成の次回の反復のための勾配として用いる。画像の再構成の各回の続く反復毎に、ブロック406において算出される実際の勾配を異なる補償項に加えて、この特定の続く反復のための勾配として用いる。一実施形態では、特定の回の反復についてブロック410で用いられる勾配は、当該特定の回の反復の全ての下位反復について勾配として用いられる。換言すると、幾つかの実施形態では、特定の回の反復についてブロック410において用いられる勾配は、当該特定の回の反復の下位反復について再計算されるのではなく、当該特定の回の反復の各回の下位反復で再利用される。
代理式という案は、任意の同時型再構成アルゴリズムに適用可能である。ここで、方法400の一実施形態を例示するOS−SPSアルゴリズムの一例について記載する。OS−SPSアルゴリズムは、データ・フィット項及び正則化項について代理式を用いる。勾配計算の計算時間を最短化するために、正則化項についての既存の代理式の冒頭にもう一つの代理式が導入され、故にこのOS−SPSアルゴリズムを「二重代理式による順序付きサブセット」アルゴリズムと呼ぶことができる。
下記の形態の一般的なPL目的関数を考える。
Φ(μ)=ΦL(μ)+ΦR(μ) (20)
式中、ΦLはデータ・フィット項であり、ΦRは正則化項である。データ・フィット項ΦLは、適当な対角行列DLと共に下記の形態の二次の代理式を有すると仮定する。
典型的な選択肢はDL=diag(dj Q,n)である。また、正則化項ΦR(μ)は、適当な対角行列DLと共に下記の形態の二次の代理式を有すると仮定する。
典型的な選択肢はDR=diag(dj R,n)である。すると、以下の「二重代理」関数が定義される。
構築によってこの二次の代理式は以下の特性を有する。すなわち、従来の正則化された順序付きサブセット法では、各々のサブセットについての最小化ステップは下式のように与えられる。
μ(n,m)=arg minμφm(μ;μ(n,m-1),μ(n,m-1)
=μ(n,m-1)
[DL(μ(n,m-1))+DR(μ(n,m-1))]-1(Dm∇ΦL(μ(n,m-1)
+∇ΦR(μ(n,m-1))) (24)
μ(n,0):=μ(n),μ(n+1):=μ(n,M) (25)
m=1,…,Mはサブセットの番号を表わす。順序付きサブセットのnによって番号付けされた各回の反復は、M回の下位反復を含む。この反復は、一つ一つのサブセット毎に正則化勾配∇ΦRを算出するため望ましくないほど遅い。この代償を少なくするために、以下の新規の更新を用いることにより二重代理式23の一般性を活用することが提案される。
μ(n,m)=arg minμφm(μ;μ(n,m-1),μ(n)
=Dφ(μ(n,m-1),μ(n))(DL(μ(n,m-1))μ(n,m-1)
+DR(μ(n))μ(n)−Dm∇ΦL(μ(n,m-1)
+∇ΦR(μ(n))) (26)
=μ(n,m-1)−Dφ(μ(n,m-1),μ(n))(Dm∇ΦL(μ(n,m-1)
+∇ΦR(μ(n)
+DR(μ(n))(μ(n,m-1)−μ(n)))…[新たな項] (27)
式中、Dφ(μ(n,m-1),μ(n))=[DL(μ(n,m-1))+DR(μ(n))]-1である。この新たな反復式27は、全てのサブセットについて同じ正則化項勾配を用いる。式24に比較して、式27の更新は正則化項勾配を更新しないことを補償する余分な項があることを除けば比較的類似している。しかしながら、この余分な項は演算が比較的単純であり(例えば加算、減算、除算、及び/又は乗算)、故に算出が速い。上の記載では、正則化項勾配は、「全ての」サブセットが更新されて初めて更新されていた。このことは、式24に記述されている本来のOS法において勾配を各回の反復毎にM回ずつ算出するのではなく、各回の反復の開始時に勾配を1回算出して、この計算を全てのサブセットについて再利用することを意味している。明らかに、正則化項勾配は必要に応じた頻度で更新されることができ、この更新頻度をUfと表わす。正則化項勾配を更新する頻度を少なくすると、早期の反復での収束速度を代償として計算費用が縮小される。
以上に述べたOS−SPSアルゴリズムを、ビュー角度を限定させた3Dコーン・ビームCT画像再構成問題について調査した。
方法400の収束速度を評価するために、n回目の反復での画像推定μ(n)と、完全収束解μ()及び真の画像μ(true)の両方との間での自乗平均平方根差を算出した。
この形式の「緩和型」OSは、最終段階で式24が収束するような1個のみのサブセットを用いるので、収束することが保証されている。図10は、164個の投影ビューによるコーン・ビームCTデータから、真のファントムの画像、FBP再構成の画像、及び完全に収束した画像(μ())を示す。
投影を41個のサブセットに分割し、すなわちサブセット当たり4個のビューに相当する。この分割は従来の選択肢に比較して寧ろ積極的な選択であって、収束を比較的顕著に高速化しようと試みる。
収束及び計算時間に対する効果を観察するために、正則化項勾配を異なる頻度で更新した。図11は、異なる正則化項更新頻度について各回の反復での収束速度を示す。図11の範囲内では、OS−41−DS−nは41個のサブセット及びUf=nによるOSを示し、「n=all」は、全てのサブセット更新が行なわれた後に1回だけ更新したことを意味する。図11(a)はμ(n)をμ()に関して示し、図11(b)はμ(n)をμ(true)に関して示している。図11は、方法400に関連するOS−SPSアルゴリズムの例が、正則化項が反復当たり比較的頻繁でなく更新されたときでも従来のOSと類似した自乗平均平方根(RMS)差を与えることを示している。
一方、図12は、同じ水準のRMS差を得るのに必要とされる計算費用が、方法400に関連するOS−SPSアルゴリズムの例によって減少したことを示している。図12は、異なる正則化項更新頻度について収束速度を計算時間に対して示している。図12(a)は、μ(n)をμ()に関して示し、図12(b)はμ(n)をμ(true)に関して示している。この例で最良の結果を与えたUf=13の場合に、方法400に関連するOS−SPSアルゴリズムの例は従来のOSよりも約3倍速く収束している。
図13は、様々な更新頻度によるOSDSの収束速度を比較している。左列において、図13は同じ時間点(初期化から4000秒後)にある画像を示している。右列において、図13は、μ()に関する絶対差画像を示している。同じ時間点での再構成画像を観察することにより、図13は、方法400に関連するOS−SPSアルゴリズムの例がより速く収束していることを明らかに示している。収束速度と及び計算費用との間にはトレードオフがあり、ここの例では、正則化項勾配を13回のサブセット更新毎に計算すると、実験に用いられた処理ハードウェアにおいては最も効率のよい結果を与えた。異なる問題については、最適な更新頻度は異なり得る。しかしながら、更新頻度を問わず、方法400に関連するOS−SPSアルゴリズムの例は従来のOSよりも速く収束することは注目に値する。問題が大規模になりサブセット数が増大するにつれて、正則化項の勾配を算出するのに必要とされる計算費用は遥かに大きくなり得る。従って、比較的実質的な利点を方法400から期待することができる。
本書に記載される方法及びアルゴリズムを用いて対象の画像を反復式で再構成する。これらの方法及びアルゴリズムは、コンピュータに記憶されて、例えば図13に示すモジュール530、ソフトウェア、ハードウェア、これらの組み合わせ、及び/又は有形の非一時的なコンピュータ可読の媒体を用いて実装される一組の命令として具現化され得る。一実施形態では、有形の非一時的なコンピュータ可読の媒体は信号を除く。
収束及び画質を改善するために係数マップを用いるという本書で用いられる概念(例えば段落[0025]で議論されたようなもの)はモデル・ベース費用関数に限らず、モデル・ベース及びモデル・ベース以外の費用関数を解く全ての同時型アルゴリズムに適用可能である。
図14は、様々な実施形態に従って形成されたイメージング・システムの一例500の見取り図である。図15は、図14に示すマルチ・モダリティ・イメージング・システム500の一部のブロック概略図である。イメージング・システムは、特に計算機式断層写真法(CT)イメージング・システム、陽電子放出断層写真法(PET)イメージング・システム、磁気共鳴イメージング(MRI)・システム、超音波イメージング・システム、X線イメージング・システム、単光子放出計算機式断層写真法(SPECT)イメージング・システム、侵襲処置用Cアーム断層写真法イメージング・システム、四肢又は乳房走査のような専用目的のCTシステム、及びこれらの組み合わせとして具現化され得る。この実施形態の例では、方法100はCTイメージング・システムに関して記載される。また、この実施形態の例では、費用関数を用いて本書に記載される様々な実施形態を説明する。
計算機式断層写真法(CT)イメージング・システムと、陽電子放出断層写真法(PET)イメージング・システム又は単光子放出計算機式断層写真法(SPECT)システムとを含む例示的な二重モダリティ・イメージング・システムの環境で様々な実施形態を説明するが、本書に記載される諸作用を果たすことが可能な他のイメージング・システムを用いることも思量されることを理解されたい。
マルチ・モダリティ・イメージング・システム500が図示されており、このイメージング・システム500はCTイメージング・システム502及びPETイメージング・システム504を含んでいる。イメージング・システム500は、異なるモダリティでの多数の走査によって単一のモダリティ・システムよりも高い診断能力を容易にすることを可能にする。一実施形態では、例示的なマルチ・モダリティ・イメージング・システム500はCT/PETイメージング・システム500である。選択随意で、CT及びPET以外のモダリティをイメージング・システム500と共に用いる。例えば、イメージング・システム500は、特に独立型CTイメージング・システム、独立型PETイメージング・システム、磁気共鳴イメージング(MRI)・システム、超音波イメージング・システム、X線イメージング・システム、及び/又は単光子放出計算機式断層写真法(SPECT)イメージング・システム、侵襲処置用Cアーム断層写真法、四肢若しくは乳房走査のような専用目的のCTシステム、並びにこれらの組み合わせであってよい。
CTイメージング・システム502はガントリ510を含んでおり、ガントリ510は、X線のビームをガントリ510の反対側に設けられた検出器アレイ514へ向けて投射するX線源512を有する。検出器アレイ514は複数の検出器素子516を含んでおり、これらの検出器素子516は横列及びチャネルを成して配列されて、被検体506のような対象を通過した投射X線を共に感知する。イメージング・システム500はまた、検出器アレイ514からの投影データを受け取って、被検体506の画像を再構成するように投影データを処理するコンピュータ520を含んでいる。動作時には、操作者が供給した命令及びパラメータをコンピュータ520によって用いて、電動式テーブル522を再配置するための制御信号及び情報を与える。さらに明確に述べると、電動式テーブル522を用いてガントリ510の内外に被検体506を移動させる。具体的には、テーブル522は、ガントリ510を通して延在するガントリ開口524を通して被検体506の少なくとも一部を移動させる。
イメージング・システム500はまた、本書に記載される様々な方法及びアルゴリズムを具現化するように構成されているモジュール530を含んでいる。モジュール530は、コンピュータ520に内蔵される単体のハードウェアとして実装され得る。選択随意で、モジュール530は、コンピュータ520にインストールされる一組の命令として実装されてもよい。この一組の命令は、独立プログラムであってもよいし、コンピュータ520にインストールされたオペレーティング・システムのサブルーチンとして組み込まれていてもよいし、コンピュータ520のインストールされたソフトウェア・パッケージの機能等であってもよい。尚、様々な実施形態は図面に示されている構成及び手段に限定されないことを理解されたい。
前述のように、検出器514は複数の検出器素子516を含んでいる。各々の検出器素子516が、被検体506を通過したときの入射X線ビームの強度を表わす、従ってビームの減弱の推定を可能にする電気信号すなわち出力を発生する。X線投影データを取得する走査中に、ガントリ510及びガントリ510に装着された構成要素は回転中心540の周りを回転する。図15は単一の横列分の検出器素子516(すなわち検出器横列1列)のみを示している。しかしながら、マルチ・スライス検出器アレイ20は、複数のスライスに対応する投影データが走査中に同時に取得され得るように、複数の平行な検出器横列を成す検出器素子516を含んでいる。
ガントリ510の回転及びX線源512の動作は、制御機構542によって制御される。制御機構542は、電力信号及びタイミング信号をX線源512に与えるX線制御器544と、ガントリ510の回転速度及び位置を制御するガントリ・モータ制御器546とを含んでいる。制御機構542に設けられたデータ取得システム(DAS)548が、検出器素子516からアナログ・データをサンプリングして、データを後の処理のためにディジタル信号へ変換する。例えば、後の処理は、本書に記載される様々な方法を具現化するためにモジュール530を用いることを含み得る。画像再構成器550が、サンプリングされてディジタル化されたX線データをDAS548から受け取って、高速画像再構成を実行する。再構成された画像はコンピュータ520に入力されて、コンピュータ520は画像を記憶装置552に記憶させる。選択随意で、コンピュータ520は、サンプリングされてディジタル化されたX線データをDAS548から受け取って、本書に記載される様々な方法をモジュール530を用いて実行することができる。コンピュータ520はまた、キーボードを有するコンソール560を介して操作者から命令及び走査パラメータを受け取る。付設されている視覚表示ユニット562によって、操作者は再構成画像及びコンピュータからの他のデータを観察することが可能になる。
操作者が供給した命令及びパラメータをコンピュータ520によって用いて、DAS548、X線制御器544及びガントリ・モータ制御器546に制御信号及び情報を与える。加えて、コンピュータ520は、電動式テーブル522を制御するテーブル・モータ制御器564を動作させて、ガントリ510に被検体506を配置する。具体的には、テーブル522は、図14に示すようにガントリ開口524を通して被検体506の少なくとも一部を移動させる。
図15へ戻り、一実施形態では、コンピュータ520は、フロッピィ・ディスク、CD−ROM、DVD、又は網若しくはインターネットのような他のディジタル・ソース等のコンピュータ可読の媒体572から命令及び/又はデータを読み取る装置570、例えばフロッピィ・ディスク・ドライブ、CD−ROMドライブ、DVDドライブ、光磁気ディスク(MOD)装置、又はイーサネット(商標)装置等のような網接続装置を含めた他の任意のディジタル装置、並びに開発中のディジタル手段を含んでいる。他の実施形態では、コンピュータ520はファームウェア(図示されていない)に記憶されている命令を実行する。コンピュータ520は本書に記載する諸作用を果たすようにプログラムされており、本書で用いられるコンピュータとの用語は当技術分野でコンピュータと呼ばれている集積回路のみに限らず、コンピュータ、プロセッサ、マイクロコントローラ、マイクロコンピュータ、プログラム可能論理コントローラ、特定応用向け集積回路、及び他のプログラム可能な回路を広範に指しており、これらの用語は本書では互換的に用いられている。
この実施形態の例では、X線源512及び検出器アレイ514は、X線ビーム574が被検体506と交差する角度が定常的に変化するように、撮像平面の内部で、撮像されている被検体506の周りをガントリ510によって回転する。一つのガントリ角度における検出器アレイ514からの一群のX線減弱測定値すなわち投影データを「ビュー」と呼ぶ。被検体506の「走査」は、X線源512及び検出器アレイ514の一回転中に異なるガントリ角度すなわちビュー角度で形成される一組のビューを含んでいる。CT走査では、投影データは、被検体506を通る二次元スライスに対応する画像を再構成するように処理される。
以上、マルチ・モダリティ・イメージング・システムの実施形態の例について詳細に説明した。図示のマルチ・モダリティ・イメージング・システムの各構成要素は本書に記載される特定の実施形態に限定されず、各々のマルチ・モダリティ・イメージング・システムの各構成要素を本書に記載される他の構成要素とは独立に別個に用いてよい。例えば、上で述べたマルチ・モダリティ・イメージング・システム構成要素を他のイメージング・システムと共に用いてもよい。
尚、様々な実施形態は、ハードウェア、ソフトウェア、又はこれらの組み合わせとして実装され得ることを特記しておく。様々な実施形態、並びに/又は構成要素、例えばモジュール、若しくは内部の構成要素及び制御器はまた、1又は複数のコンピュータ又はプロセッサの一部として具現化され得る。コンピュータ又はプロセッサは、計算装置、入力装置、表示ユニット、及びインタフェイス例えばインターネットにアクセスするためのインタフェイスを含み得る。コンピュータ又はプロセッサはマイクロプロセッサを含み得る。マイクロプロセッサは通信バスに接続され得る。コンピュータ又はプロセッサはまた、メモリを含み得る。メモリは、ランダム・アクセス・メモリ(RAM)及び読み出し専用メモリ(ROM)を含み得る。コンピュータ又はプロセッサはさらに、記憶装置を含んでいてよく、記憶装置はハード・ディスク・ドライブ、又は固体ドライブ及び/若しくは光ドライブ等のような着脱自在の記憶ドライブであってよい。記憶装置はまた、コンピュータ又はプロセッサにコンピュータ・プログラム又は他の命令を読み込む他の同様の手段であってよい。
本書で用いられる「コンピュータ」との用語は、マイクロコントローラ、縮小命令セット・コンピュータ(RISC)、特定応用向け集積回路(ASIC)、論理回路、GPU、FPGA、及び本書に記載された作用を果たすことが可能な他の任意の回路又はプロセッサを用いたシステムを含む任意のプロセッサ方式又はマイクロプロセッサ方式のシステムを含み得る。上の例は例示のみのためのものであり、従って「コンピュータ」との語の定義及び/又は意味を限定しないものとする。コンピュータ又はプロセッサは、入力データを処理するために1又は複数の記憶要素に記憶されている一組の命令を実行する。記憶要素はまた、所望又は必要に応じてデータ又は他の情報ソースを記憶していてよい。記憶要素は、処理機械の内部の情報ソース又は物理メモリ素子の形態にあってよい。
上述の一組の命令は、本発明の様々な実施形態の方法及び工程のような特定の動作を実行するように処理機械であるコンピュータ又はプロセッサに命令する様々な命令を含み得る。一組の命令は、ソフトウェア・プログラムの形態にあってよい。ソフトウェアは、システム・ソフトウェア又はアプリケーション・ソフトウェアのような様々な形態にあってよく、非一時的なコンピュータ可読の媒体であってよい。さらに、ソフトウェアは、別個のプログラムの集合、より大きなプログラムの内部のプログラム・モジュール又はプログラム・モジュールの一部の形態にあってよい。ソフトウェアはまた、オブジェクト指向プログラミングの形態のモジュール型プログラミングを含み得る。コンピュータ又はプロセッサによる入力データ処理は、利用者の命令に応答して行なわれてもよいし、以前の処理の結果に応答して行なわれてもよいし、他の処理機械によって発行された要求に応答して行なわれてもよい。
本書で用いる場合には、単数形で記載されており単数不定冠詞を冠した要素又はステップとの用語は、排除を明記していない限りかかる要素又はステップを複数備えることを排除しないものと理解されたい。さらに、本発明の「一実施形態」に対する参照は、所載の特徴を同様に組み入れている追加の実施形態の存在を排除しないものと解釈されたい。また、反対に明記されていない限り、特定の特性を有する1若しくは複数の要素を「含んでいるcomprising」又は「有しているhaving」実施形態は、この特性を有しないような追加の要素も包含し得る。
また、本書で用いられる「画像を再構成する」との表現は、画像を表わすデータが生成されるが可視画像は形成されないような本発明の実施形態を排除するものではない。従って、本書で用いられる「画像」との用語は可視画像及び可視画像を表わすデータの両方を広く指す。但し、多くの実施形態は少なくとも1枚の可視画像を形成し又は形成するように構成されている。
本書で用いられる「ソフトウェア」及び「ファームウェア」との用語は互換的であり、RAMメモリ、ROMメモリ、EPROMメモリ、EEPROMメモリ、及び不揮発性RAM(NVRAM)メモリを含めたメモリに記憶されて、コンピュータによって実行される任意のコンピュータ・プログラムを含んでいる。以上のメモリ形式は例示のみのためのものであり、従ってコンピュータ・プログラムの記憶に利用可能なメモリの形式に関して制限するものではない。
以上の記載は例示説明のためのものであって制限するものではないことを理解されたい。例えば、上述の各実施形態(及び/又は各実施形態の諸観点)を互いに組み合わせて用いてよい。加えて、本発明の範囲を逸脱することなく、特定の状況又は材料を発明の教示に合わせて適応構成する多くの改変を施すことができる。本書に記載されている材料の寸法及び形式は、発明の各パラメータを定義するためのものであるが、限定するものではなく例示する実施形態である。以上の記載を吟味すれば、当業者には他の多くの実施形態が明らかとなろう。従って、本発明の範囲は、特許請求の範囲に関連して、かかる特許請求の範囲が網羅する等価物の全範囲と共に決定されるものとする。特許請求の範囲では、「including包含する」との用語は「comprising含む」の標準英語の同義語として、また「in whichこのとき」との用語は「whereinここで」の標準英語の同義語として用いられている。また、特許請求の範囲では、「第一」、「第二」及び「第三」等の用語は単にラベルとして用いられており、これらの用語の目的語に対して数値的要件を課すものではない。さらに、特許請求の範囲の制限は、「手段プラス機能(means-plus-function)」形式で記載されている訳ではなく、かかる特許請求の範囲の制限が、「〜のための手段」に続けて他の構造を含まない機能の言明を従えた文言を明示的に用いていない限り、合衆国法典第35巻第112条第6パラグラフに基づいて解釈されるべきではない。
この書面の記載は、最適な態様を含めて発明の様々な実施形態を開示し、また任意の装置又はシステムを製造して利用すること及び任意の組み込まれた方法を実行することを含めてあらゆる当業者が発明の様々な実施形態を実施することを可能にするように実例を用いている。特許付与可能な発明の様々な実施形態の範囲は特許請求の範囲によって画定されており、当業者に想到される他の実例を含み得る。かかる他の実例は、特許請求の範囲の書字言語に相違しない構造要素を有する場合、又は特許請求の範囲の書字言語と非実質的な相違を有する等価な構造要素を含む場合には、特許請求の範囲内にあるものとする。
100、200、300、400:対象の画像を反復式で再構成する方法
500:イメージング・システム
502:CTイメージング・システム
504:PETイメージング・システム
506:被検体
510:ガントリ
512:X線源
514:検出器アレイ
516:検出器素子
520:コンピュータ
522:電動式テーブル
524:ガントリ開口
530:モジュール
540:回転中心
542:制御機構
544:X線制御器
546:ガントリ・モータ制御器
548:データ取得システム(DAS)
550:画像再構成器
552:記憶装置
560:コンソール
562:視覚表示ユニット
564:テーブル・モータ制御器
570:命令及び/又はデータを読み取る装置
572:非一時的なコンピュータ可読の媒体
574:X線ビーム

Claims (28)

  1. 対象の画像を反復式で再構成する方法であって、
    前記画像に関連する測定データを入手するステップと、
    前記画像を再構成するために同時型アルゴリズムを用いるステップと
    を備えており、前記画像を再構成するために前記同時型アルゴリズムを用いる前記ステップは、
    ボクセル依存型であるスケーリング・ファクタを決定するステップと、
    前記画像を再構成するために前記ボクセル依存型スケーリング・ファクタを目的関数の勾配に適用するステップと
    を含んでいる、方法。
  2. ボクセル依存型であるスケーリング・ファクタを決定する前記ステップは、取得幾何学的構成に基づいて直接前記スケーリング・ファクタを算出することを含んでいる、請求項1に記載の方法。
  3. ボクセル依存型であるスケーリング・ファクタを決定する前記ステップは、取得幾何学的構成に基づいて直接前記スケーリング・ファクタを算出することを含んでおり、前記測定データは三次元(3D)取得幾何学的構成において非一様にサンプリングされている、請求項1に記載の方法。
  4. ボクセル依存型であるスケーリング・ファクタを決定する前記ステップは、係数マップを生成するステップと、前記ボクセル依存型スケーリング・ファクタを算出するために前記係数マップを用いるステップとを含んでいる、請求項1に記載の方法。
  5. 係数マップを生成する前記ステップは、少なくとも2回の反復又は少なくとも2回の下位反復の一方についての差から得られる更新を変換することを含んでいる、請求項4に記載の方法。
  6. 係数マップを生成する前記ステップは、局所的に滑らかであり低域通過フィルタ処理を施されることが可能な係数マップを生成することを含んでいる、請求項4に記載の方法。
  7. 係数マップを生成する前記ステップは、初期画像から前記係数マップを得ることを含んでいる、請求項4に記載の方法。
  8. 係数マップを生成する前記ステップは、初期画像に関する事前情報、現在の画像の特徴、又は取得幾何学的構成の少なくとも一つにより前記係数マップを強化することを含んでいる、請求項4に記載の方法。
  9. 前記画像を再構成するために前記ボクセル依存型スケーリング・ファクタを目的関数の勾配に適用する前記ステップは、前記ボクセル依存型スケーリング・ファクタを用いて前記勾配をスケーリングすることを含んでいる、請求項1に記載の方法。
  10. 前記スケーリングは、ATWAを用いて得られる又は一つ一つのボクセルへのサブセットの寄与に関係付けられるの少なくとも一方である、請求項9に記載の方法。
  11. 前記画像を再構成するために前記ボクセル依存型スケーリング・ファクタを目的関数の勾配に適用する前記ステップは、前記測定データの各々のスライス毎に緩和係数を算出することを含んでおり、前記測定データは三次元(3D)取得幾何学的構成において非一様にサンプリングされており、着目領域(ROI)の内部のスライスは、前記ROIの外部のスライスに比較して異なる数のビューによりサンプリングされている、請求項1に記載の方法。
  12. 前記画像を再構成するために前記ボクセル依存型スケーリング・ファクタを目的関数の勾配に適用する前記ステップは、各回の反復毎ではなくn回の反復毎まで減少させた頻度で前記ボクセル依存型スケーリング・ファクタを更新することを含んでいる、請求項1に記載の方法。
  13. 前記係数マップはダウンサンプリングされている、請求項4に記載の方法。
  14. 検出器アレイと、
    該検出器アレイに結合されたコンピュータと
    を備えたイメージング・システムであって、前記コンピュータは、
    対象の画像に関連する測定データを取得し、
    前記画像を再構成するために同時型アルゴリズムを用いる
    ように構成されており、前記コンピュータは、ボクセル依存型であるスケーリング・ファクタを決定して、前記画像を再構成するために前記ボクセル依存型スケーリング・ファクタを目的関数の勾配に適用することにより、前記画像を再構成するために前記同時型アルゴリズムを用いるように構成されている、イメージング・システム。
  15. 前記コンピュータは、取得幾何学的構成に基づいて直接前記スケーリング・ファクタを算出することにより、ボクセル依存型である前記スケーリング・ファクタを決定するように構成されている、請求項14に記載のイメージング・システム。
  16. 前記コンピュータは、係数マップを生成して、前記ボクセル依存型スケーリング・ファクタを算出するために前記係数マップを用いることにより、ボクセル依存型である前記スケーリング・ファクタを決定するように構成されている、請求項14に記載のイメージング・システム。
  17. 前記コンピュータは、少なくとも2回の反復又は少なくとも2回の下位反復の一方についての差から得られる更新を変換することにより、前記係数マップを生成するように構成されており、前記変換は離散型又は連続型の一方である、請求項16に記載のイメージング・システム。
  18. 前記コンピュータは、前記ボクセル依存型スケーリング・ファクタを用いて前記勾配をスケーリングすることにより、前記ボクセル依存型スケーリング・ファクタを前記目的関数の前記勾配に適用するように構成されている、請求項14に記載のイメージング・システム。
  19. 前記コンピュータは、前記測定データの各々のスライス毎に緩和係数を算出することにより、前記ボクセル依存型スケーリング・ファクタを前記目的関数の前記勾配に適用するように構成されており、前記測定データは三次元(3D)取得幾何学的構成において非一様にサンプリングされており、着目領域(ROI)の内部のスライスは、前記ROIの外部のスライスに比較して異なる数のビューによりサンプリングされている、請求項14に記載のイメージング・システム。
  20. 対象の画像を反復式で再構成する方法であって、
    前記画像に関連する測定データを入手するステップと、
    前記画像を再構成するために同時型アルゴリズムを用いるステップと
    を備えており、前記同時型アルゴリズムは正則化関数の勾配を算出することを含んでおり、前記画像を再構成するために前記同時型アルゴリズムを用いる前記ステップは、
    複数のサンプル勾配値をルックアップ・テーブルに記憶させるステップと、
    前記ルックアップ・テーブルの前記サンプル勾配値を入手するステップと、
    前記サンプル勾配値を用いて近似的な実際の勾配値を算出するステップと
    を含んでいる、方法。
  21. 前記サンプル勾配値を用いて近似的な実際の勾配値を算出する前記ステップは、前記再構成の各回の反復毎に近似的な実際の勾配値を算出することを含んでいる、請求項20に記載の方法。
  22. 前記サンプル勾配値を用いて近似的な実際の勾配値を算出する前記ステップは、前記近似的な実際の勾配値を算出するために補間方式を用いることを含んでいる、請求項20に記載の方法。
  23. 複数のサンプル勾配値をルックアップ・テーブルに記憶させる前記ステップは、固定されたサンプリング周波数で前記勾配値をサンプリングすることを含んでいる、請求項20に記載の方法。
  24. 対象の画像を反復式で再構成する方法であって、
    前記画像に関連する測定データを入手するステップと、
    前記画像を再構成するために同時型アルゴリズムを用いるステップと
    を備えており、前記同時型アルゴリズムは正則化関数の勾配を算出することを含んでおり、前記画像を再構成するために前記同時型アルゴリズムを用いる前記ステップは、
    前記画像の前記再構成の反復において実際の勾配を算出するステップと、
    前記実際の勾配に補償項を加えるステップと、
    前記画像の前記再構成の次回の反復のための前記勾配として前記実際の勾配及び前記補償項を用いるステップと
    を含んでいる、方法。
  25. 前記画像の前記再構成の次回の反復のための前記勾配として前記実際の勾配及び前記補償項を用いる前記ステップは、前記次回の反復に続く前記画像の前記再構成の各回の反復毎に異なる補償項を前記実際の勾配に加えることを含んでいる、請求項24に記載の方法。
  26. 前記画像の前記再構成の次回の反復のための前記勾配として前記実際の勾配及び前記補償項を用いる前記ステップは、前記次回の反復の全ての下位反復について前記勾配を用いることを含んでいる、請求項24に記載の方法。
  27. 前記実際の勾配に補償項を加える前記ステップは、前記正則化関数についての代理式を構築することを含んでいる、請求項24に記載の方法。
  28. 前記補償項は加算、減算、除算、又は乗算の少なくとも一つを含むのみである、請求項24に記載の方法。
JP2013127082A 2012-06-22 2013-06-18 反復式再構成の方法及び装置 Active JP6212294B2 (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US13/531,255 US8958660B2 (en) 2012-06-22 2012-06-22 Method and apparatus for iterative reconstruction
US13/531,255 2012-06-22

Publications (3)

Publication Number Publication Date
JP2014004359A true JP2014004359A (ja) 2014-01-16
JP2014004359A5 JP2014004359A5 (ja) 2016-08-04
JP6212294B2 JP6212294B2 (ja) 2017-10-11

Family

ID=49486632

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2013127082A Active JP6212294B2 (ja) 2012-06-22 2013-06-18 反復式再構成の方法及び装置

Country Status (5)

Country Link
US (1) US8958660B2 (ja)
JP (1) JP6212294B2 (ja)
CN (1) CN103514615B (ja)
DE (1) DE102013106466A1 (ja)
NL (1) NL2011023C2 (ja)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017006640A (ja) * 2015-06-18 2017-01-12 株式会社東芝 X線コンピュータ断層撮影装置、逐次近似再構成方法および医用画像処理装置
WO2017154217A1 (ja) * 2016-03-11 2017-09-14 株式会社島津製作所 画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置
JP2017189612A (ja) * 2016-04-11 2017-10-19 東芝メディカルシステムズ株式会社 放射線画像診断装置及び医用画像処理装置
JP2018044957A (ja) * 2016-09-13 2018-03-22 キヤノンメディカルシステムズ株式会社 医用画像処理装置、医用画像診断装置及び医用画像処理方法
JP2018046976A (ja) * 2016-09-21 2018-03-29 株式会社島津製作所 逐次近似画像再構成方法、逐次近似画像再構成プログラムおよび断層撮影装置
JP2019535009A (ja) * 2016-09-30 2019-12-05 コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. ノイズに起因するアーチファクト形成の動的抑制を用いた反復画像再構成
JP2020521134A (ja) * 2017-05-22 2020-07-16 プリズマティック、センサーズ、アクチボラグPrismatic Sensors Ab 画像再構築のための方法および装置
CN112489149A (zh) * 2019-08-23 2021-03-12 西门子医疗有限公司 用于重建医学图像数据的计算机实现的方法

Families Citing this family (24)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9508164B2 (en) * 2013-01-23 2016-11-29 Czestochowa University Of Technology Fast iterative image reconstruction method for 3D computed tomography
WO2014160766A1 (en) * 2013-03-26 2014-10-02 The Johns Hopkins University Task-based source-detector trajectories for tomographic imaging
US9495770B2 (en) 2013-08-14 2016-11-15 University Of Utah Research Foundation Practical model based CT construction
EP3123446B1 (en) * 2014-03-26 2020-07-08 Koninklijke Philips N.V. Image generation apparatus
EP3123399A4 (en) 2014-03-27 2017-10-04 Hrl Laboratories, Llc System for filtering, segmenting and recognizing objects in unconstrained environments
CN104331293A (zh) * 2014-11-07 2015-02-04 江苏理工学院 一种基于fpga的n阶数据拟合方法及n阶数据拟合器
CN104637033B (zh) * 2014-12-30 2018-02-06 深圳先进技术研究院 Ct内部感兴趣区域成像方法和系统
CN107592935A (zh) * 2015-03-20 2018-01-16 伦斯勒理工学院 X射线ct的自动系统校准方法
US9721361B2 (en) * 2015-05-29 2017-08-01 General Electric Company Systems and methods for parallel processing of imaging information
DE102015219622A1 (de) * 2015-10-09 2017-04-13 Siemens Healthcare Gmbh Rekonstruktion einer Abbildung anhand einer oder mehrerer Bildgebungsmodalitäten
EP3365872B1 (en) * 2015-10-20 2019-12-11 Koninklijke Philips N.V. Device and method for reconstructing an x-ray computed tomography image
CN108292428B (zh) * 2015-12-11 2023-03-31 上海联影医疗科技股份有限公司 图像重建的系统和方法
KR102539407B1 (ko) * 2016-01-14 2023-06-01 프리스매틱 센서즈 에이비 x-선 검출기용 측정 회로, 및 해당 방법 및 x-선 영상 시스템(A MEASUREMENT CIRCUIT FOR AN X-RAY DETECTOR, AND A CORRESPONDING METHOD AND X-RAY IMAGING SYSTEM)
EP3245634B1 (en) * 2016-03-31 2020-10-14 Shanghai United Imaging Healthcare Co., Ltd. System and method for image reconstruction
US10332252B2 (en) * 2016-12-29 2019-06-25 General Electric Company Slope constrained cubic interpolation
US10692251B2 (en) * 2017-01-13 2020-06-23 Canon Medical Systems Corporation Efficient variance-reduced method and apparatus for model-based iterative CT image reconstruction
TWI639414B (zh) * 2017-11-17 2018-11-01 財團法人國家同步輻射研究中心 電腦斷層掃描影像的對位方法
US20190180481A1 (en) * 2017-12-13 2019-06-13 General Electric Company Tomographic reconstruction with weights
CN108460241B (zh) * 2018-05-28 2022-08-02 中国民用航空中南地区空中交通管理局 一种仪表着陆系统扰动仿真方法
CN110470743B (zh) * 2019-08-23 2021-11-16 天津大学 电学/超声信息融合的双模态层析成像方法
CN111932648B (zh) * 2020-06-17 2023-05-12 北京信息科技大学 一种由螺旋采样光场数据重建三维物体的方法
CN112244884B (zh) * 2020-10-27 2023-08-29 沈阳先进医疗设备技术孵化中心有限公司 骨图像获取方法、装置、控制台设备及ct系统
US20220261964A1 (en) * 2021-02-12 2022-08-18 Flir Commercial Systems, Inc. Image non-uniformity mitigation systems and methods
CN116580115B (zh) * 2023-03-29 2023-12-05 赛诺威盛医疗科技(扬州)有限公司 螺旋ct图像迭代重建方法、装置及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS61225640A (ja) * 1985-03-30 1986-10-07 Toshiba Corp 断層撮影装置
US20030156684A1 (en) * 2002-02-20 2003-08-21 Fessler Jeffrey A. Method for statistically reconstructing images from a plurality of transmission measurements having energy diversity and image reconstructor apparatus utilizing the method
US20100177946A1 (en) * 2007-05-18 2010-07-15 Marleen De Bruijne Semi-automatic contour detection
US20110142314A1 (en) * 2009-12-15 2011-06-16 Jiang Hsieh System and method of increasing temporal resolution of an x-ray image
US20110164031A1 (en) * 2010-01-06 2011-07-07 Kabushiki Kaisha Toshiba Novel implementation of total variation (tv) minimization iterative reconstruction algorithm suitable for parallel computation

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006116316A2 (en) * 2005-04-22 2006-11-02 University Of Chicago Open source trajectory method and apparatus for interior imaging
US7760848B2 (en) * 2006-09-08 2010-07-20 General Electric Company Method and system for generating a multi-spectral image of an object
JP2010527741A (ja) * 2007-05-31 2010-08-19 ジェネラル エレクトリック カンパニー 画像再構成において利得変動の補正を容易にする方法及びシステム
US9171353B2 (en) * 2008-07-16 2015-10-27 Siemens Medical Solutions Usa, Inc. Multimodal image reconstruction with zonal smoothing
US8600136B2 (en) * 2008-09-19 2013-12-03 Koninklijke Philips N.V. Method for generation of attenuation map in PET-MR
US8761478B2 (en) * 2009-12-15 2014-06-24 General Electric Company System and method for tomographic data acquisition and image reconstruction
US8189735B2 (en) * 2010-07-22 2012-05-29 General Electric Company System and method for reconstruction of X-ray images

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS61225640A (ja) * 1985-03-30 1986-10-07 Toshiba Corp 断層撮影装置
US20030156684A1 (en) * 2002-02-20 2003-08-21 Fessler Jeffrey A. Method for statistically reconstructing images from a plurality of transmission measurements having energy diversity and image reconstructor apparatus utilizing the method
US20100177946A1 (en) * 2007-05-18 2010-07-15 Marleen De Bruijne Semi-automatic contour detection
JP2010527647A (ja) * 2007-05-18 2010-08-19 ノルディック・バイオサイエンス・イメージング・アクティーゼルスカブ 半自動式輪郭検出方法
US20110142314A1 (en) * 2009-12-15 2011-06-16 Jiang Hsieh System and method of increasing temporal resolution of an x-ray image
JP2011125700A (ja) * 2009-12-15 2011-06-30 General Electric Co <Ge> X線画像の時間分解能を高めるシステム及び方法
US20110164031A1 (en) * 2010-01-06 2011-07-07 Kabushiki Kaisha Toshiba Novel implementation of total variation (tv) minimization iterative reconstruction algorithm suitable for parallel computation
JP2011139894A (ja) * 2010-01-06 2011-07-21 Toshiba Corp 画像処理方法及びx線コンピュータ断層撮影装置

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017006640A (ja) * 2015-06-18 2017-01-12 株式会社東芝 X線コンピュータ断層撮影装置、逐次近似再構成方法および医用画像処理装置
WO2017154217A1 (ja) * 2016-03-11 2017-09-14 株式会社島津製作所 画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置
JPWO2017154217A1 (ja) * 2016-03-11 2019-01-10 株式会社島津製作所 画像再構成処理方法、画像再構成処理プログラム並びにそれを搭載した断層撮影装置
JP2017189612A (ja) * 2016-04-11 2017-10-19 東芝メディカルシステムズ株式会社 放射線画像診断装置及び医用画像処理装置
JP7055606B2 (ja) 2016-09-13 2022-04-18 キヤノンメディカルシステムズ株式会社 医用画像処理装置、医用画像診断装置及び医用画像処理方法
JP2018044957A (ja) * 2016-09-13 2018-03-22 キヤノンメディカルシステムズ株式会社 医用画像処理装置、医用画像診断装置及び医用画像処理方法
JP2018046976A (ja) * 2016-09-21 2018-03-29 株式会社島津製作所 逐次近似画像再構成方法、逐次近似画像再構成プログラムおよび断層撮影装置
JP2019535009A (ja) * 2016-09-30 2019-12-05 コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. ノイズに起因するアーチファクト形成の動的抑制を用いた反復画像再構成
JP7041128B2 (ja) 2016-09-30 2022-03-23 コーニンクレッカ フィリップス エヌ ヴェ ノイズに起因するアーチファクト形成の動的抑制を用いた反復画像再構成
JP7041128B6 (ja) 2016-09-30 2022-05-30 コーニンクレッカ フィリップス エヌ ヴェ ノイズに起因するアーチファクト形成の動的抑制を用いた反復画像再構成
JP2020521134A (ja) * 2017-05-22 2020-07-16 プリズマティック、センサーズ、アクチボラグPrismatic Sensors Ab 画像再構築のための方法および装置
JP7143330B2 (ja) 2017-05-22 2022-09-28 プリズマティック、センサーズ、アクチボラグ 画像再構築のための方法および装置
CN112489149A (zh) * 2019-08-23 2021-03-12 西门子医疗有限公司 用于重建医学图像数据的计算机实现的方法

Also Published As

Publication number Publication date
CN103514615B (zh) 2018-01-30
JP6212294B2 (ja) 2017-10-11
US20130343673A1 (en) 2013-12-26
CN103514615A (zh) 2014-01-15
DE102013106466A1 (de) 2014-05-22
NL2011023C2 (en) 2014-10-07
NL2011023A (en) 2013-12-24
US8958660B2 (en) 2015-02-17

Similar Documents

Publication Publication Date Title
JP6212294B2 (ja) 反復式再構成の方法及び装置
JP6280700B2 (ja) 反復式再構成の方法、非一時的なコンピュータ可読の媒体及びイメージング・システム
JP5530637B2 (ja) 画像再構成の方法及びシステム
US7885371B2 (en) Method and system for image reconstruction
US9020230B2 (en) Method and apparatus for iterative reconstruction
Ritschl et al. Improved total variation-based CT image reconstruction applied to clinical data
Hauptmann et al. Multi-scale learned iterative reconstruction
JP5543448B2 (ja) 高効率コンピュータ断層撮影
US9489752B2 (en) Ordered subsets with momentum for X-ray CT image reconstruction
JP6275826B2 (ja) ノイズ除去再構成画像データエッジ改善
US9600924B2 (en) Iterative reconstruction of image data in CT
JP6169558B2 (ja) コントラスト依存の解像度をもつ画像
US10692251B2 (en) Efficient variance-reduced method and apparatus for model-based iterative CT image reconstruction
WO2015189811A1 (en) Improved image reconstruction for a volume based on projection data sets
Jin et al. A model-based 3D multi-slice helical CT reconstruction algorithm for transportation security application
JP2018020120A (ja) 医用画像処理装置及び医用画像処理プログラム
CN113424227A (zh) 锥束计算机断层扫描(cbct)中的运动估计和补偿
JP7362460B2 (ja) 医用画像処理装置、方法及び記憶媒体
US10515467B2 (en) Image reconstruction system, method, and computer program
US10307114B1 (en) Iterative volume image reconstruction using synthetic projection images
WO2018022565A1 (en) System and method for tomographic image reconstruction
Brendel et al. Penalty weighting for statistical iterative CT reconstruction

Legal Events

Date Code Title Description
A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20160614

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20160614

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20170220

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20170228

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20170525

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

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20170915

R150 Certificate of patent or registration of utility model

Ref document number: 6212294

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