JP5902686B2 - パラレル磁気共鳴撮像法ならびにダイナミックおよびパラレル磁気共鳴撮像法の実施方法 - Google Patents

Description

この発明は、パラレルの、ファンクショナルＭＲＩ（ｆＭＲＩ）などの動的（時間分解）磁気共鳴法を含む、人体のパラレル磁気共鳴法（ｐＭＲＩ）の実施方法に関する。

全体の収集時間（global acquisition time）を減少させることは、医用磁気共鳴法（ＭＲＩ）における一番の関心事であり、さらにｆＭＲＩのような動的撮像法に関するときはなおさらである。実際には、短い収集時間が、収集したｆＭＲＩデータの空間／時間分解能を改善することを可能にし、より効率的な統計的分析をもたらす。加えて、全体の撮像時間を減少させることによって、患者の動きによって生じるいくつかのさらなるアーチファクト（artifact）を回避することができる。このような理由から、パラレル撮像法システムが開発されている：周波数領域（すなわち、いわゆるｋ空間）中で、少なくとも１つの空間方向、すなわち位相変調方向に沿ったナイキストサンプリングレートよりもＲ倍低いレートでサンプリングされたデータを同時に収集するために、下に置かれた被写体または人体の周囲に位置する補完的な感度特性を有する複数受信表面コイルが用いられる；ここで、Ｒは通常“減少係数（reduction factor）”と呼ばれる。それゆえ、全体の収集時間は、従来型の非パラレル撮像法によるものよりもＲ倍短い。その後、個々の受信部によって得られたアンダーサンプリング（undersampled）された“基本（elementary）”画像を展開することによって完全な視野（ＦＯＶ）を構築するために再構成手順が実施される。アンダーサンプリングレートに関連した偽信号のアーチファクトに起因するパラレルＭＲＩ（ｐＭＲＩ）内の低信号対雑音比（ＳＮ比）のため、再構成は困難だがやりがいのある課題である。これらの影響は、収集プロセス中のノイズに起因し、コイル感度マップの推定中の誤差の存在にも起因する。

空間高調波の同時収集（ＳＭＡＳＨ）［Sodicksonら1997］は、ｋ空間領域で動作する最初の再構成法であった。それは欠落した位相変調手順を生成する、事前推定されたコイル感度マップの線形結合を用いる。

再構成法に基づくいくつかの他のｋ空間は、GRAPPA（一般化自己較正部分パラレル収集）［Griswoldら2002］、そしてＳＥＮＳＥ（感度符号化）［Pruessmannら1999］のようにも提案されている。ＳＥＮＳＥは、第１に縮小されたＦＯＶ像の再構成に依存し、第２に空間展開技術に依存する２段階法であり、加重最小二乗法に相当する。この技術は、リファレンススキャン（通常２Ｄグラジェントエコー（ＧＲＥ））を用いたコイル感度マップの正確な推定を要する。それが現在最も頻繁に用いられ、特に頭部および心臓撮像法に応用されるｐＭＲＩ技術である。
ｐＭＲＩ中の再構成法の総括には［Hogeら2005］を参照せよ。

ＳＥＮＳＥはしばしばノイズのない完全なコイル感度マップ情報の場合において、正確な再構成を実現するものと考えられており、上述の全ての方法にも当てはまる。しかしながら、実際問題として、データ中のノイズの存在とコイル感度マップの推定の不正確さが、再構成問題を悪い状態にする。

画像再構成は不良逆問題であり、正則化法は一般的に完全なＦＯＶ像をより良く推定するために用いられる。これらの技術の大部分は、画像領域内で動作する；特に、これはTikhonov正則化［Yingら2004］の場合であり、これは、滑らかな制約を促進し、また再構成画像と事前のリファレンス画像の間の差の二乗を説明する２次のペナルティ（penalty）項を用いる。しかしながら、（1.5テスラまでの）低い磁場強度が用いられるとき、正則化の使用にも関わらず、（R=2の値を超える）高い減少係数が通常、実現不可能であると考えられる。なぜなら、再構成画像が深刻な偽信号のアーチファクトを受けるためである。

［Chaariら2008］と［Chaariら2009］において、この発明者は、ウェーブレット基底正則化方式を用いたパラレルMRI中の正則化画像再構成を実施する方法を記載しており、減少係数Rの増大を可能にする。

この発明は、前記方法のいくつかの改善を提供することを目的とし、それを動的撮像法（例えばｆＭＲＩ）にまで拡張し、完全または部分的に自己較正（または“管理なし（unsupervised）”）にする。

それゆえ、この発明の目的は、
− 既知のまたは推定された感度マップおよび雑音共分散行列を有するそれぞれの受信アンテナから前記人体の一連の基本磁気共鳴画像であって、ｋ空間でアンダーサンプリングされた基本画像を収集する工程と、
− 前記人体の磁気共鳴画像の正則化再構成を行う工程とを備え、
前記磁気共鳴画像の正則化再構成を行う工程は、
− 収集した前記基本画像に基づき再構成された再構成画像の尤度を表す誤差項および
− 前記再構成画像のフレーム係数の実際の統計的分布と前記係数の事前分布との間の偏差を表すフレームペナルティ項
を備えた評価関数を最小化することによって離散フレーム空間で実施され、
再構成画像のフレーム係数の前記事前分布は、前記人体の補助的な磁気共鳴画像に基づいて推定される、人体のパラレル磁気共鳴撮像法である。

前記磁気共鳴撮像法の別の実施態様によれば、
− 前記誤差項は、前記収集基本画像に基づき再構成された前記再構成画像の負の対数尤度を表すことができ、
− 前記人体の磁気共鳴画像の正則化再構成を行う前記工程は、前記収集基本磁気共鳴画像とフレーム係数の前記事前分布を考慮して、人体の画像を規定する一連のフレーム係数の完全な事後分布を、前記フレーム空間において、最大化することによって行うことができ、
− 前記人体の前記補助的な磁気共鳴画像は、前記収集基本磁気共鳴画像から再構成でき、

特に、前記人体の前記補助的な磁気共鳴画像は、感度エンコーディング−SENSE−アルゴリズムを用いて前記収集基本磁気共鳴画像から再構成でき、さらに、前記人体の前記補助的な磁気共鳴画像は、非正則化SENSEアルゴリズム；画像空間で正則化されたSENSEアルゴリズム；およびk空間において正則化されたSENSEアルゴリズム：から選択されたアルゴリズムを用いた前記収集基本磁気共鳴画像から再構成でき、
− 前記フレーム係数の汎用ガウス−ラプラス事前統計的分布を想定し、そして前記分布のパラメータは、最尤推定量または事後平均推定量を用いて、前記人体の前記補助的な磁気共鳴画像に基づき推定でき、
− 前記誤差項は、二次平均誤差項となり得、
− 前記収集基本画像は三次元像となり得、そして磁気共鳴画像の正則化再構成を行う前記工程は、離散三次元フレーム空間で実施され得、
− さらに、前記収集三次元基本画像は、画像を収集すべき被験者をスライスした二次元基本画像を積層させることによって得ることができ、磁気共鳴画像の− 正則化再構成を行う前記工程は、冗長ウェーブレットフレーム表現に基づくことができ、
− 或いは、磁気共鳴の正則化再構成を行う前記工程は、非冗長なウェーブレット表現に基づくことができ、
− 前記評価関数は、再構成画像の全変分ノルム；および凸制約から選択された少なくとも１つの空間領域ペナルティ項も備えることができ、

この発明の他の目的は、
− 既知または推定された感度マップおよび雑音共分散行列を有するそれぞれの受信アンテナから前記人体の一連の時系列の基本磁気共鳴画像であって、ｋ空間でアンダーサンプリングされた基本画像を収集する工程と、
− 前記人体の一時系列の磁気共鳴画像の正則化再構成を行う工程とを備え、
一時系列の基本磁気共鳴画像の正則化再構成を行う工程は、
− 対応する収集基本画像に基づき再構成された各再構成画像の尤度を表す誤差項；および
− その系列の連続した画像間の画素ごとまたはボクセルごとの差を表す、時間的ペナルティ項
を備えた評価関数を最小化することによって実施される、人体のダイナミックおよびパラレル磁気共鳴撮像法の実施方法である。

前記方法の別の実施態様によれば：
− 前記時間的ペナルティ項が、エッジ保存関数、特に凸エッジ保存関数、そしてさらにｐ≧１、好ましくは１≦ｐ<１．５のＬ_pノルムに基づくことができ、
前記時間的ペナルティ項は、第１部分時間的ペナルティ項および第２部分時間的ペナルティ項の和によって与えられ：第１部分時間的ペナルティ項は、その系列の各偶数画像と奇数画像との間のピクセルごとまたはボクセルごとの差を表し；そして、第２の部分時間的ペナルティ項は、その系列の各奇数画像と先行する偶数画像との間のピクセルごとまたはボクセルごとの差を表し；前記評価関数は、前記第１および第２の部分時間的ペナルティ項に対し、近接写像を用いることによって最小化される。
− 磁気共鳴画像の正則化再構成を行う前記工程が、離散フレーム空間で実施でき、前記評価関数は、前記係数の各再構成画像のフレーム係数の統計的分布と前記係数の事前分布との間の偏差を示す、フレームペナルティ項をも備えることができ；再構成画像のフレーム係数の前記事前分布は、前記人体の補助的な磁気共鳴画像に基づき推定される。この場合において、前記基本画像は三次元画像となり得、前記離散フレーム空間は、離散三次元フレーム空間となり得る。
− その方法は、最尤推定量を用いて；前記時間的ペナルティ項の重み付けパラメータを自動的に決定する工程をさらに有することができ、さらに正確には、その方法は、再構成すべき画像の、各ピクセル若しくはボクセルまたは隣接したピクセル若しくはボクセルに対する時間的ペナルティ項の前記重み付けパラメータを推定することができ、さらに、前記時間的ペナルティ項がL_pノルムに基づくことができ、時間的ペナルティ項の前記重み付けパラメータおよびパラメータpが前記最尤推定量を用いて共に決定できる。特に、時間的ペナルティ項の前記重み付けパラメータおよびパラメータｐは、制約ｐ≧１の下で前記最尤推定量を用いて共に決定できる。
− 前記誤差項は、基準とみなされる基本磁気共鳴画像に対する前記基本磁気共鳴画像のそれぞれの剛体変換を規定する幾何学的パラメータに依存でき、そして一時系列の基本磁気共鳴画像の正則化再構成を行う前記工程は、前記幾何学的パラメータについて前記関数を最小化することによっても実施される。
− 前記基本画像は、エコープラナー撮像法によって収集できる。
− 前記基本画像は、４より大きいかまたは４に等しい減少係数でアンダーサンプリングできる。

この発明のもう１つの特徴と有利な点は、添付図面に関連した次の説明から明らかになるだろう：

図１は、人間の脳の磁気共鳴画像のウェーブレット係数の実部および虚部の実験に基づいたヒストグラムである。 図２は、R=2で、先行技術により知られたＴｉｋｈｏｎｏｖ正則化SENSE pMRI再構成法を用いて、得られた人間の脳の９枚の解剖学的な軸位断スライスである。 図３は、R=4で、先行技術により知られたＴｉｋｈｏｎｏｖ正則化SENSE pMRI再構成法を用いて、得られた人間の脳の９枚の解剖学的な軸位断スライスである。 図４は、R=2（左パネル）およびR=4（右パネル）で、先行技術により知られたTV-正則化SENSE pMRI再構成法を用いて、得られた同一の人間の脳の解剖学的な軸位断スライスである。 図５は、R=2で、この発明の実施形態による自動較正２Ｄウェーブレット変換に基づく正則化方式（UWR-SENSE）を用いて、得られた同一の人間の脳の９枚の解剖学的な軸位断スライスである。 図６は、R=4で、この発明の実施形態による自動較正２Ｄウェーブレット変換に基づく正則化方式（UWR-SENSE）を用いて、得られた同一の人間の脳の９枚の解剖学的な軸位断スライスである。 図７は、R=2で、この発明の別の実施形態による自動較正２Ｄウェーブレット変換に基づく正則化方式（CWR-SENSE）を用いて得られた同一の人間の脳の９枚の解剖学的な軸位断スライスである。 図８は、R=4で、この発明の別の実施形態による自動較正２Ｄウェーブレット変換に基づく正則化方式（CWR-SENSE）を用いて得られた同一の人間の脳の９枚の解剖学的な軸位断スライスである。 図９は、R=4で、この発明の別の実施形態による自動較正結合ウェーブレット全変分正則化方式を用い、そして非冗長の正規直交ウェーブレット基底上への画像分解を用いて、得られた同一の人間の脳の９枚の解剖学的な軸位断スライスである。 図１０は、R=4で、この発明の別の実施形態による自動較正結合ウェーブレット全変分正則化方式を用い、そして２つの正規直交基底の結合によって構成された冗長ウェーブレットフレーム上への画像分解を用いて、得られた同一の人間の脳の９枚の解剖学的な軸位断スライスである。 図１１は、従来の非パラレルＭＲＩ（上段）、この発明の別の実施形態による、R=4でmSENSEパラレルＭＲＩ再構成（中段）および自動較正３Ｄウェーブレット変換に基づく正則化方式（３Ｄ-UWRSENSE）（下段）：を用いて得られた、同一の人間の脳の９枚の解剖学的な軸位断スライスである； 図１２は、従来の非パラレルＭＲＩ（上段）、この発明の別の実施形態による、R=4でmSENSEパラレルＭＲＩ再構成（中段）および自動較正３Ｄウェーブレット変換に基づく正則化方式（３Ｄ-UWR-SENSE）（下段）：を用いて得られた、同一の人間の脳の３枚の解剖学的な矢状断スライスである； 図１３は、TV正則化（左列）；非冗長ウェーブレット基底上への画像分解を用いたこの発明の実施形態による結合ウェーブレット全変分正則化方式（中列）および２つの正規直交基底の結合によって構成した冗長ウェーブレットフレーム上への画像分解を用いたこの発明の実施形態による結合ウェーブレット全変分正則化方式（右列）：を用いて得られた同一の人間の脳の３枚の解剖学的な軸位断スライス（上段）およびそれらの拡大詳細図（下段）である； 図１４は、エコープラナー（EPI）ｆＭＲＩを用いて収集し、先行技術により知られた、mSENSE方式およびこの発明の実施形態による自動較正２Ｄウェーブレット変換に基づく正則化方式（４Ｄ-UWR-SENSE）を用いて再構成した、人間の脳の軸位断、冠状断および矢状断スライスである； 図１５は、ＥＰＩ ｆＭＲＩを用いて検出されたaC-aS対比の、解剖学的なイメージングに重ね合わせた、被験者レベルのスチューデント−ｔマップ；データは、それぞれR = 2（図の上部）およびR = 4（図の下部）で、mSENSE、UWR-SENSEおよび４Ｄ-UWR-SENSE方式を用いて再構成した；矢状断、冠状断および軸位断像を示す； 図１６は、R = 2で、それぞれmSENSE、UWR-SENSEおよび４Ｄ-UWR-SENSE方式を用いて再構成した２人の被験者に対するaC-aS対比の、解剖学的なイメージングに重ね合わせた、被験者レベルのスチューデント−ｔマップである；矢状断、冠状断および軸位断像を示す； 図１７は、EPI ｆＭＲＩを用いて検出したLc-Rc対比の、解剖学的なイメージングに重ね合わせた、被験者レベルのスチューデント−ｔマップ；データは、それぞれR = 2（図の上部）およびR = 4（図の下部）で、mSENSE、UWR-SENSEおよび４Ｄ-UWR-SENSE方式を用いて再構成した；矢状断、冠状断および軸位断の図を示す； 図１８は、R = 4で、それぞれmSENSE、UWR-SENSEおよび４Ｄ-UWR-SENSE方式を用いて再構成した２人の被験者に対するLc-Rc対比の被験者レベルのスチューデント−ｔマップである；矢状断、冠状断および軸位断像を示す； 図１９は、R = 2およびR=4について、mSENSE、UWR-SENSEおよび４Ｄ-UWR-SENSEを用いてデータを再構成した、aC-aS対比についてのグループレベルのスチューデント−ｔマップである；矢状断、冠状断および軸位断像を示す； 図２０は、R = 2およびR=4について、mSENSE、UWR-SENSEおよび４Ｄ-UWR-SENSEを用いてデータを再構成した、Lc-Rc対比についてのグループレベルのスチューデント−ｔマップである；矢状断、冠状断および軸位断像を示す；

この発明の詳細を述べる前に、ｐＭＲＩ（そして、特にＳＥＮＳＥ）について、いくつかの基本的事実を思い起こす必要があるだろう。
ＭＲＩは、２Ｄ（２次元）または直接３Ｄ（３次元）のいずれでも実施することができる撮像技術であり、複雑なＲＦパルス設計に依存する。２次元の場合、隣接するスライスを用いて体積を取り扱う。ＰＭＲＩは、その限定が（すなわち２または３次元内の）どのようなものであっても、この方法がk空間の走査をより速くすることから、いずれの状況にも適用し得る。簡単のため、２Ｄの場合を以下に示す。

調査中の被写体（例えば、患者の頭の脳）内のスピン密度ρを測定するため、Ｃ個のコイルの配列が用いられる。画像収集は特定の撮像手順に基づく。この発明の実施形態において、解剖学的ＭＲＩが３Ｄ ＭＰＲＡＧＥ手順に基づく一方で、ファンクショナルＭＲＩは２Ｄエコープラナー撮像法（ＥＰＩ）を用いて実施される。以下、２Ｄの場合に焦点を合わせて説明する；各コイルl（１≦ｃ≦Ｃ）が受信する信号ｄ_cは、k空間内のいくつかの位置ｋ＝（ｋ_y，ｋ_x）^T（上付記号^Tは、転置を意味する）において評価されたコイルの感度プロフィールｓ_cによって加重した所望の２Ｄ場ρのフーリエ変換である。
受信信号、

はそれゆえ以下のサンプリング方式によって特徴付けられる：

は、付加白色ガウシアン雑音（ＡＷＧＮ）を具体化させたものであり、ｒ＝（ｙ，ｘ）^Tは画像領域内の空間位置である。簡単のため、直交座標系が通常採用される。最も簡単な形式で、ＳＥＮＳＥイメージングがフーリエ変換の分離可能性に起因する１次元逆問題を解決する。Δｙ＝Ｙ／Ｒをサンプリング周期とする。ここで、Ｙは位相変調方向に沿った視界（ＦＯＶ）のサイズであり、ｙを同一方向に沿った画像領域内の位置とし、ｘを周波数エンコード方向に沿った画像領域内の位置とし、Ｒ≦Ｌを減少係数（reduction factor）とする。２Ｄ逆フーリエ変換は、空間領域内の測定信号の復元を可能にする。k空間のアンダーサンプリングをＲによって説明することによって、式（１）を以下の行列形式に表現し直すことができる。

式（２）において、（ｎ（ｒ））_rは、回転ゼロ平均ガウシアン複素値ランダムベクトルの列である。これらの雑音ベクトルは、i.i.d.（独立の同一分布）であり、サイズＣ×Ｃの共分散行列Ψに関して空間的に独立である。Ψは、高周波パルスを有しない全てのコイルから得られるデータによって推定される。そして、２つのコイルｃ₁およびｃ₂の間の共分散に相当する一般式Ψ（ｃ₁，ｃ₂）は、次式によって与えられる：

ここで、（・）^*は、複素共役を表す。統計モデルとしては、コイル間の統計的独立または最近接統計的独立モデルのような行列Ψを想定できる。最初の場合において、行列Ψは対角行列となり、コイル従属変数

は、［Ｄｏｎｏｈｏ，１９９５］によって記述された平均絶対偏差（ＭＡＤ）技術を用いて推定することができ、これは付加白色ガウシアン雑音（ＡＷＧＮ）という前提に依存する。最も簡単な場合においては、Ψは恒等行列に等しいとみなす。ゼロでない非対角項が存在するとき、式（４）によって与えられる実験に基づく共分散が、これらの非対角項を推定するために用いられる。

従って、再構成手順は、逆変換（２）と空間位置ｒ＝（ｙ，ｘ）^Tにおけるｄ（ｒ）からρバー（ｒ）を復元することからなる。｜ρバー｜は、視覚化のためにのみ考慮されるが、データ（ｄ_c）_1≦l≦Cおよび未知の画像ρは複素値であることに注意する。

ＳＥＮＳＥ法［Ｐｒｕｅｓｓｍａｎｎ他、１９９９］とも呼ばれる簡単な再構成法は、加重最小二乗（ＷＬＳ）を基準とした最小化に基づくものである。その目的は、各空間位置ｒにおけるベクトルρハット_WLS（ｒ）を見つけることである：

であり、（・）^Hは転置複素共役を表し、（・）^#は疑似逆を表し、そして以下の式、

は、Ｃ^L空間上のノルムを定義する。

上で議論したように、この逆問題は、一般に不良設定問題であり、正則化、例えば、Ｔｉｋｈｏｎｏｖ正則化を要する。この正則化過程は、通常ρハット_PWLS（ｒ）を計算して、以下のペナルティ重み付き最小二乗（ＰＷＬＳ）を基準として最小化することにある：

ここで、Ｉ_RはＲ次元恒等行列である。その正則化パラメータκ＞０は、データに近い値とペナルティ項との間のバランスを確保する。ペナルティ項は所定の基準ベクトルρ_r（ｒ）からのずれを調整する。解ρハット_PWLS（ｒ）は、以下の閉形式の式を与える：

この解の正確さは、基準ベクトルおよび正則化パラメータκの選択に依存することに留意する。

Ｔｉｋｈｏｎｏｖ正則化は、画像中のボケ（blurring）を発生させることが知られている。
このような限界を克服すべく、画像領域内に適用され、ボケ効果を制限することによって正則化をより効率化して画像境界を保存する“エッジ保存”ペナルティ項が提案されている。しかしながら、これらの項の導入は、常に数値的に解くことが容易でない微分不可能な最適化問題を生じることがある。

上述のように、［Chaariら2008］および［Chaariら2009］において、この発明者はウェーブレット変換（WT）に基づく；新たな正則化方式を提案しており、この明細書において要約される；
実際には、SENSEに基づき再構成された画像において、十分に空間的に局在した不自然な結果が非常に高い、または低い輝度で歪んだ曲線として現れ、WTは有用な情報の良好な空間および周波数の局所化を可能にする強力な道具として認識されている。
ウェーブレットおよびウェーブレット変換の概略紹介が［Pesquet-Popescu, Pesquet］によって提供されている。

以下において、TはWT作用素を表す。それは、jmaxの分解能のレベルを超えて実施される可分の２Ｄ Mバンドウェーブレット基底上への離散分解に対応する。
サイズY×Xの対象画像ρバーは、標準的な内積＜・｜・＞およびノルム||・||を有するユークリッド空間、C^K空間の要素とみなすことができる。(e_k)_1≦k≦Kを、C^K空間の離散ウェーブレット基底とみなす。ウェーブレット分解作用素Ｔは、線形作用素として定義される：

次に、再構成目的に役立つ共役作用素Ｔ^*は、全単斜線形作用素として定義される。

対象画像関数の結果として生じるウェーブレット係数体は、ζ = ((ζ_a,k)_{1≦k≦Kjmax}, (ζ_0,j,k)_{1≦k≦Kjmax,1≦k≦Kj})によって定義される。ここで、K_j = KM^-2j は（YおよびXがM^jmaxの倍数であると想定して）分解能ｊで所定のサブバンド中のウェーブレット係数の数である。そして係数は、ζ_a,kが分解能レベルj_maxでの近似係数を示し、ζ_0,j,kが分解能レベルjおよび方向o∈O=｛0，…，M - 1｝²＼｛(0,0)｝での詳細係数を示すようにインデックスが再構成されている。
２進の場合（Ｍ＝２）、水平、垂直または斜め方向に対応する３方向がある。正規直交ウェーブレット基底を考慮したとき、共役作用素Ｔ^*は、逆WT作用素T^-1に変換し、Tの作用素ノルム||T||は１に等しい。
対象画像ρの推定は、再構成ウェーブレット作用素Ｔ^*によって行われる。ζバーをρバー=T^*ζバーのような未知のウェーブレット係数とする。目的は、観測dから係数ζバーのベクトルの推定量ζハットを作ることである。これは、ウェーブレット係数の適切な事前分布に依存するベイズ的アプローチに基づく。
式（２）において観測モデルと雑音に関する条件（ゼロ平均i.i.d.ガウス分布コイル間の相関行列Ψ）を仮定すると、尤度関数はYX視野内のピクセル上で因数分解できる。

ここで、p = T^*ζバーであり、

ｆをウェーブレット領域中の像の事前確率密度関数（pdf）とする。ここで、ウェーブレット係数の実部および虚部が独立であると仮定する。ウェーブレット係数の実（虚）部は、各サブバンドにおいて、互いに独立で同一の分布に従うことも仮定する。しかしながら、これらの統計的特性は、２つの異なるサブバンド間で変わることがある。さらに、考慮されたウェーブレット係数の実部および虚部の経験分布を見ることによって、それらの実験に基づいたヒストグラムが“一般化ガウス−ラプラス”（GGL）分布によく合うことに、この発明者は気付いている。GGLは、単一モードを示し、その形状がガウス密度とラプラス密度との間で変化する。対応するpdfは、以下のように記載される：

ここでα∈R₊空間とβ∈R_-空間は、推定すべきハイパーパラメータである。図１は、長さ８のDaubechiesフィルタで２進（Ｍ＝２）ウェーブレット分解を用いた一次分解能レベルで水平方向の詳細なサブバンドの実部と虚部の実験に基づいたヒストグラムを例示する。この図は、採用されたGGL分布が一般化ガウス（GG）pdfより実験に基づいたヒストグラムにより良く適合することも示す。

最も粗い分解能レベルj_maxにおいて、近似係数の実部と虚部のいずれの分布も低周波数のサブバンドに属するため、ガウス分布によってモデル化できる。

種々の選択が可能だが、その親しみやすさと簡易さから、MAP（最大事後確率）推定量が推定目的で用いられる。上で与えられた事前分布および尤度に基づき、MAP推定量が完全事後分布の最大化または負の対数尤度の最小化によって計算され、

または以下の基準を最小化することによって同等に計算される。

ここで、Re(・)およびIm(・)(または・^Reおよび・^Im)は、それぞれ実部および虚部を表す。事前パラメータ（ハイパーパラメータ）α_0,j=(α_0,j ^Re, α_0,j ^Im)∈(Ｒ^* ₊空間)², β_0,j=(β_0,j ^Re, β_0,j ^Im)∈(Ｒ^* ₊空間)², μ_0,j=(μ_0,j ^Re, μ_0,j ^Im)∈(Ｒ空間)² およびσ_0,j=(σ_0,j ^Re, σ_0,j ^Im)∈(Ｒ₊空間)² は、未知であり推定する必要がある。ハイパーパラメータの推定は、この発明の重要な部分であり、後で広範囲にわたって議論する。

J_Lはリプシッツ連続勾配で微分可能である一方で、J_pは微分可能でない。それゆえ、ペナルティ関数J_WT は凸関数であるが、その最適化手法は疑似共役勾配のような従来の凸最適化技術に依存し得ない。

[Daubechiesら2004; Chauxら2007]において開発された反復最適化手法の一般形を用いて最適化を行うことができる。
関数、

ここで、φ^Reとφ^1mは、Γ₀(R^K空間)内の関数であり、Re(x)(Im(x))は、近接演算子が次のように定義される、x∈C^K空間の要素の実部（虚部）のベクトルである。

[Chauxら2007]におけるアルゴリズムを複数の場合に拡張することによって、Ｊ_WTの極小値が以下のアルゴリズム（アルゴリズム１）に従って反復的に計算し得る。ここで、Ｊ_Lの勾配を最初に計算し、次にフレーム係数を更新する。λ_nとγ_nがそれぞれ緩和およびステップサイズパラメータに対応することに注目し得る。

アルゴリズム１：
（γ_n）_n>0と（λ_n）_n>0は、正の実数列である。
１：ｎ＝０およびε≧０とおく。ζ⁽ⁿ⁾を初期化してＪ⁽ⁿ⁾＝Ｊ_WT（ζ⁽ⁿ⁾）とおく。
２：repeat
３：ρ⁽ⁿ⁾＝Ｔ^*ζ⁽ⁿ⁾とおくことによって画像を再構成する。
４：以下のように画像ｕ⁽ⁿ⁾を計算する：
∀ｒ∈｛１，…，Ｙ／Ｒ｝×｛１，…，Ｘ｝，
ｕ⁽ⁿ⁾（ｒ）＝２Ｓ^H（ｒ）Ψ^-1（Ｓ（ｒ）ρ⁽ⁿ⁾（ｒ）−ｄ（ｒ）），
ここで、ベクトルｕ⁽ⁿ⁾（ｒ）は、ρバー（ｒ）がρバーから定義されるのと同様にｕ⁽ⁿ⁾から定義される（（２．１７）式を参照）。
５：ｕ（ｎ）のウェーブレット係数ｖ（ｎ）＝Ｔｕ⁽ⁿ⁾＝（ｖａ，（ｖ_0,j）_{o∈O空間,1≦j≦jmax}）を決定する。
６：再構成された画像ρ⁽ⁿ⁺¹⁾の近似係数を更新する：
∀ｋ∈｛１，…，Ｋ_jmax｝，ζ_a,k ⁽ⁿ⁺¹⁾＝ζ_a,k ⁽ⁿ⁾＋λ_n（ｐｒｏｘ_γnΦa（ζ_a,k ⁽ⁿ⁾−γ_nｖ_a,k ⁽ⁿ⁾）−ζ_a,k ⁽ⁿ⁾）
７：再構成された画像ρ⁽ⁿ⁺¹⁾の詳細係数を更新する：
∀ｏ∈Ｏ空間，∀ｊ∈｛１，…，ｊ_max｝，∀ｋ∈｛１，…，Ｋ_j｝，ζ_o,j,k ⁽ⁿ⁺¹⁾＝ζ_o,j,k ⁽ⁿ⁾＋λ_n（ｐｒｏｘ_γnΦo,j（ζ_o,j,k ⁽ⁿ⁾−γ_nｖ_o,j,k ⁽ⁿ⁾）−ζ_o,j,k ⁽ⁿ⁾）
８：Ｊ⁽ⁿ⁺¹⁾＝Ｊ_WT（ζ⁽ⁿ⁺¹⁾）
９：ｎ←ｎ＋１
１０：until ｜Ｊ⁽ⁿ⁾−Ｊ^(n-1)｜≦εＪ^(n-1)
１１：return ρ⁽ⁿ⁾＝Ｔ^*ζ⁽ⁿ⁾

各ｒ∈｛１，…，Ｙ｝×｛１，…，Ｘ｝につき、θ_r≧０をエルミート正半定符号行列

の最大固有値とし、θ＝ｍａｘ_{r∈｛1,…,Y｝×｛1,…,X｝}θ_r＞０とする。アルゴリズム１の収束を担保すべく、ステップサイズおよび緩和パラメータが以下の条件に従う：

[Chaariら2009]に議論されるように、最適化問題(11), (12)は付加（好ましくは凸の）制約、例えば、画像強度値の上界および下界を含むように修正される。
この発明は、上記のウェーブレット正則化を用いることにより、再構成のアーチファクトを平滑化する一方で、画像の詳細を提供することが可能になるが、画像の均質な領域中にいくつかの不規則を生じることがある。一方、他の（非ウェーブレット基底の）正則化方式は、平滑領域を正則化するために適用されることが知られているが、画像の詳細の過度の平滑化を犠牲にする。これは特に、とりわけ[Rajら2007]および[Blockら2007]によって記載された“全変分”ＴＶ正則化の場合である。
この発明の一面を構成する、[Chaariら2008]および[Chaariら2009]の方法の第１の改善は、WTとTVの異なる特性を活用するため、結合正則化の枠組みでWTとTVを結合して、それらに相互の欠点を軽減させることにある。結合ウェーブレット変換の全変分（WTTV）正則化は、以下のペナルティ付き基準の最適化に帰着できる。

ここで、Ｊ_LおよびＪ_Pは上で定義されたものであり、κ₁＞０およびκ₂＞０は、正則化パラメータであり、Ｔ^*はウェーブレット随伴作用素である。画像ρ∈C空間^Y×Xの離散全変分ノルムは、次式によって与えられる：

ここで、各ρ∈C空間^Y×Xに対し、∇₁は、水平方向に平滑化された勾配演算子であり、次式によって定義される。

説明のため、画像ρが周期的または画像の境界がトロイダル形状と同等であると仮定する。
２項より多くの項が関係するため、最適性基準(15)を最小化することは、(11)-(13)によって与えられる空間のみの基準を最小化するよりもずっと困難である。その一方で、上記のフォワード−バックワード(FB)アルゴリズムだけが、２項からなる基準の最小化に適用される。２項より多くの項からなる微分不可能な凸関数の反復最小化は、いわゆる並行近接アルゴリズム（PPXA）[CombettesおよびPesquet,2008]を用いて実施でき、３つの関連する項のそれぞれの近接演算子の計算をも要する。
問題は、以下の近接演算子の計算から生じる。

このような問題を回避するため、式(15)のTVペナリゼーションを[Combettes and Pesquet, 2008]のように４項に分割できる。このTVペナルティ項はそれゆえ以下のように記載される。

ここで、各(q,r)∈[0,1]²に対し、

ここで、[0,1]中のqおよびrのそれぞれに対し、↓_q,rは、以下の式によって定義されるデシメーション（decimation）演算子とする。

そして、Ｕ_q+2rを以下の演算子とする。

ここで、各ρ∈C空間^Y×Xに対し、∇₀および∇₂は、次式によって定義される。

また、ｈを次式によってC空間^Y×Xに定義される関数であるとする。

また、上式から以下の式が分かる。

その結果、(15)式の最適化問題は、以下のように書くことができる：

正則化WT-TV再構成法は、アルゴリズム2に要約される。ここで、PPXAアルゴリズム[Combettes and Pesquet, 2008]がJ=6の凸関数からなる式(18)の最適性基準を最小化するために用いられる。

アルゴリズム2

となるように、γ∈]0, + ∞ [, n = 0, (ω_i)_1≦i≦6 ∈ [0,1]⁶を設定する。ここで、全てのi∈｛1,…,6｝に対し、

である。また、ε≧0と定め、

およびJ⁽ⁿ⁾=0と初期化する。
１：repeat
２：

となるように画像u⁽ⁿ⁾を計算する。ここで、ρ₁ ⁽ⁿ⁾=T^*ζ₁ ⁽ⁿ⁾である。
３：u⁽ⁿ⁾のウェーブレット係数p₁ ⁽ⁿ⁾=Tu⁽ⁿ⁾を計算する
４：

を計算する
５：全てのi∈｛0,1,2,3｝に対し、

を計算する
６：

とおく
７：λ_n∈[0,2]と設定する
８：for i=1 to 6 do
９：以下

１０：end for
１１：ζ⁽ⁿ⁺¹⁾=ζ⁽ⁿ⁾+λ_n(P⁽ⁿ⁾-ζ⁽ⁿ⁾)
１２：J⁽ⁿ⁺¹⁾=J_WT-TV(ζ⁽ⁿ⁺¹⁾)を計算する
１３：ｎ←ｎ＋１
１４：until |J⁽ⁿ⁾-J^(n-1)|≦εJ^(n-1)
１５：return ρ⁽ⁿ⁾=T^*ζ⁽ⁿ⁾の値を返す

[Chaariら2008]および[Chaariら2009]の方法において、３次元（３Ｄ）画像化は、正則化された２次元（２Ｄ）画像を積み重ねることによって実施される。この発明のもう１つの側面を構成する、この方法のさらなる改善点は、直接の、正則化された３Ｄ画像の再構成にある。３Ｄの場合において、画像の再構成の問題はさらに、以下のように書かれる。

これは、式(2)と同様である；しかしながら、r = (y, x, z)は現在3次元の空間位置であり、z∈[1, ... ,Z]は第３の方向に沿った位置（スライス選択）である。さらに、ペナルティ項Ｊ_pは、３Ｄの２進ウェーブレット変換のウェーブレット係数の分布に依存する。

この発明のもう１つの側面を構成する、この方法のさらなる改善点は、４次元(４Ｄ)の、動的ＭＲＩの時空正則化の実施にある。この発明の実施形態は、例えば脳のｆＭＲＩに適用することができ、４Ｄのデータセットを得るには、脳全体の体積を数回収集しなければならない。従来のｆＭＲＩにおいて、３Ｄ画像は同一のｆＭＲＩセッションに属するが、これらは独立であると想定される。しかしながら、実際には、３Ｄ時間画像はどういうわけか独立である。なぜなら、それらは同一の実験的パラダイムに作用する同一のｆＭＲＩセッションに属するからである。
ＢＯＬＤ（血中酸素濃度依存の）時系列および収集雑音は、実際には脳のｆＭＲＩセッションにおいて時間的な相関がある。このような理由のため、３Ｄ画像間の時間依存性を考慮に入れることは、得られた画像全体のＳＮ比の増加に役立ち、それゆえ、ｆＭＲＩの統計的分析の信頼性を高める。しかしながら、動的ＭＲＩにおいて、画像物体形状は一般に収集中に変化することから、時間的正則化に再構成過程を結びつけることは、非常に困難である。

N_r個の３次元画像（N_rは通常、偶数）の時系列の４Ｄ再構成を処理するため、(12)式および(19)式の観測モデルが以下のように書かれる：

ここで、t∈[1, . . . , N_r]は、収集時間である。２進３Ｄウェーブレット作用素Ｔを用いて、係数が以下のように再構成される

追加の時間的c_p正則化項を考慮して、４Ｄ“体積”（すなわち、３Ｄ画像または体積の時系列）の再構成が以下の最適性基準の最小化によって実施される。

ここで、ζ = (ζ¹，ζ² ，. . . ，ζ^Nr )T、κ>0は正則化パラメータ、J_pは(13)式のように定義される。
(21)式において、

は、対応する基本画像が得られた場合、それぞれ再構成された画像の尤度を示す“誤差項”である。

は、列の連続した画像間のピクセルごと、またはボクセルごとの差を示す、時間的ペナルティ項である。

は、空間的正則化のために用いられるウェーブレットペナルティ項である。異なる（例えば、非ウェーブレット基底の）空間的正則化項も用いられることが、理解されよう。

他の収集画像に関する時間依存性を考慮することによって、再構成された３Ｄ画像ρ^tを収集時間tで得るため、随伴ウェーブレット作用素T^*がその後ζの各要素ζ^tに適用される。

時間的ペナルティ項は、ウェーブレットを基底としないが、列の連続した画像間のピクセルごと、またはボクセルごとの差を示す。さらに正確には、(21)式において、時間的ペナルティ項は、Ｌ_pノルムに基づく；この発明者は、pが1より大きく、または1に等しいこと、そして好ましくは、1≦p<1.5であることを見出している。特に、列の連続した画像間のピクセルごと、またはボクセルごとの差の（好ましくは、凸の）エッジ保存関数を基底とする、時間的ペナルティ項の他の形式が用いられ得る。

最適性基準(21)の最小化は、(11)-(13)による空間のみの基準を最小化するよりもずっと難しい。なぜなら、２項以上の項が関わる一方で、上記のフォワード−バックワード（FB）アルゴリズムが２項からなる基準の最小化にのみ適用されるからである。２項以上の項を含む微分不可能な凸関数の反復最小化は、既に引用されたパラレルプロクシマルアルゴリズム（PPXA）[CombettesおよびPesquet,2008]を用いて実施でき、３つの関連する項のそれぞれの近接演算子の計算も要する。この作業は、(21)式の２つの第１項に対しては、かなり単純である。なぜなら、それらの変数は、時間の変数tおよび空間位置rに関して分離可能だからである。しかしながら、これは、対応する近接演算子の非自明な計算をする時間ペナルティ項（(21)式の第３項）の場合には当てはまらない。
その後、時間的ペナルティを時間変数ｔ（所定の収集時間ｔは第１項または第２項のいずれにも含まれる）に関して分離可能であり、近接写像の計算が容易な２項に分解することによって、J_STの最適性基準を書き換えることが提案されている。より具体的には、時間的ペナルティ項J_Tは、第１の部分時間的ペナルティ項J_T ¹と第２の部分時間的ペナルティ項J_T ²との和として表現され、ここで：
第１の部分時間的ペナルティ項は、列の偶数番目の画像と先行する奇数番目の画像との間のピクセルごと、またはボクセルごとの差を表すものである。
そして第２の部分時間的ペナルティ項は、列の奇数番目の画像と先行する偶数番目の画像との間のピクセルごと、またはボクセルごとの差を表すものである。

J_T ¹およびJ_T ²は、時間変数ｔに関して分離可能であり、対応する近接演算子は関係する項の各々の近接演算子に基づき容易に計算できる。以下の関数を検討する。

そして、Ｈは次式によって定義される線形作用素である。

関連する随伴作用素Ｈ^*はそれゆえ次式によって与えられる。

Φの近接写像は、次式によって与えられる：

その結果得られる空間−時間的最適性基準J_STの最小化のためのアルゴリズム（アルゴリズム３）は、以下で与えられる。

アルゴリズム３

となるように、γ∈]0, +∞[, n = 0, (ω_i)_1≦i≦4 ∈ [0,1]⁴を定める。ここで、全てのi∈｛1,…,4｝およびt∈｛1,…,N_r｝に対し、

である。また、ε≧0と定め、

およびJ⁽ⁿ⁾=0と初期化する。
１：repeat
２：p₄ ^1,(n)=ζ₄ ^1,(n).
３：for t=1 to N_r do
４：

となるように画像u_t ⁽ⁿ⁾を計算する。ここで、ρ_t,1 ⁽ⁿ⁾=T^*ζ₁ ^t,(n)である。
５：ウェーブレット係数p₁ ^t,(n)=Tu_t ⁽ⁿ⁾を計算する
６：

を計算する
７： if tが偶数 then
８：

を計算する。
９： else if tが奇数かつt>1 then
１０：

を計算する。
１１： end if
１２： if t>1 then
１３：

とおく。
１４： end if
１５： end for
１６： p₄ ^Nr,(n)=ζ₄ ^Nr,(n)とおく。
１７：

とおく。
１８：

とおく。
１９：λ_n ∈ [0,2]と設定する
２０：for i=1 to 4
２１：以下

２２：end for
２３：ζ⁽ⁿ⁺¹⁾=ζ⁽ⁿ⁾+λ_n(P⁽ⁿ⁾-ζ⁽ⁿ⁾)
２４：J⁽ⁿ⁺¹⁾=J_ST(ζ⁽ⁿ⁺¹⁾)を計算する
２５：ｎ←ｎ＋１
２６：until |J⁽ⁿ⁾-J^(n-1)|≦εJ^(n-1)
２７：ζハット = ζ⁽ⁿ⁾とおく。
２８：return t∈｛１，…，N_r｝に対し、ρハット_t = T^*ζハット^t

空間−時間正則化は、たとえあまり一般的でないとしても、２次元画像の時系列の場合にも適用できる。さらに、空間的ペナルティ項はまた、純空間的な２Ｄまたは３Ｄ正則化の場合におけるように、全変分成分を含むこともある。

４ＤのｆＭＲＩデータセットにおいて、それは注目に値し、撮像される体の体積を何度も収集するため、前記体（すなわち、脳）は連続するスキャンの間にわずかに動くことがある。その固有の動きアーチファクトは、ｆＭＲＩ分析にとって劇的であることもある。このような状況において、そして標準的なｆＭＲＩ研究が厳格な変形を適用することによる運動補正手順を含むことから、推定法の他の拡張が提案されている。

この方法は、推定段階における剛体変形パラメータδｒを取り込むことによって、結果として生ずる動きアーチファクトを説明する。この変形は、推定段階のために用いられる標準のSENSEデータに最初に適用される上記方法の拡張も提案されており、標準のSENSE再構成データではないが、平滑マトリックス（すなわち、QR-SENSE）を有するTikhonov正則化を用いた２次的に正則化された（QR）版のSENSEに関する推定手順を実施する。QR-SENSEは、SENSEにおける追加のフィルタリング手順を導入するより、推定手順のロバスト性を改善するため、より効率的であるように思われ、閉形式の式を認めることから計算効率も良いままである。ｆＭＲＩにおいて発生する動きアーチファクトをさらに考慮するため、提案された再構成手法（アルゴリズム３）が、再構成過程の間、剛体変化パラメータを説明するために拡張することもできる。式(22)中の最適性基準は、それゆえ以下のように書き換えられる：

ここで、運動補正の転換を適用し、以下を想定することによって、ζ_c ^tがζ^tから導かれる

J_T ¹（J_T ²）は、時間変数ｔに関して分離可能であり、その近接演算子が上記総和中の関係項のそれぞれの近接演算子に基づき容易に計算できる。

剛体変形パラメータδｒは、これらのパラメータについて22-2式によって表されるJ_ST(ζ)を最小化することによって見出すことができる。代わりに、前記剛体変形パラメータは、従来のSENSE再構成に基づき、別々に推定することもできる。

ここに記載された全ての正則化法（ウェーブレット基底２Ｄ空間正則化、ウェーブレット基底３Ｄ空間正則化、ハイブリッドウェーブレット基底／全変分２Ｄおよび３Ｄ空間正則化、３Ｄおよび４Ｄ空間−時間正則化）は、適切な推定アルゴリズムを用いて手動または自動のいずれでも設定されるべき数値パラメータを含む。これらのパラメータは、ウェーブレット係数の空間分布の先行するパラメータ、すなわち“ハイパーパラメータ”を含む。

それらは、時間的正則化パラメータκ、いわゆる“形状パラメータ”、すなわち、時間的正則化（21式参照）で用いられるLpノルム||・||_pを定義する“ｐ”の値、およびウェーブレットの相対的重さおよびハイブリッド正則化における全変分ペナルティ項の相対的重みを決定するパラメータκ₁、κ₂を含むこともできる。

既に引用された文献[Chaariら2008] および [Chaariら2009]は、相対パラメータをどのように決定するかを教示しない。この発明の他の面を構成するこれらの文献によって記載された方法のさらなる改善点は、自動較正アルゴリズムを提供することにあり、完全または部分的な非管理ウェーブレット基底正則化を可能にする。分解のタイプに依存するこの目標を実現すべく、２つの異なる枠組みが考慮される。正規直交ウェーブレット分解の場合、ハイパーパラメータは最尤法（以下を参照）を用いて正確に推定できる。対照的に、冗長フレーム分解をとるとき、確率的サンプリング手順がより効率的である。後者の場合において、ハイパーパラメータは、最小平均二乗誤差（MMSE）または同等である事後平均推定量を用いて推定される（以下、参照）。事後平均推定量は、これがフレームの特別な場合であるとき、正規直交ウェーブレット分解で用いることもできる。

第１に、２Ｄまたは３Ｄ空間的正則化のために用いられるウェーブレット係数の事前確率分布を定義するハイパーパラメータの決定が考慮される。

ハイパーパラメータを決定するための第１の確率

は、以下の式に従い、積分された尤度を最大化するものである。

Cが未知のため、（24）式を最大化することは、欠測データ問題である。
求められた画像分解ζを完全に統合し、[Dempsterら1977]によって記載された集中EMアルゴリズムを用いて画像再構成およびハイパーパラメータ推定間で反復することが必要である。計算の負荷を軽減するため、参照または“補助”全視界画像ρバーが利用可能であることを仮定することによって、異なって続行させることが有利であり、ウェーブレット分解ζバー=(T^*)^-1ρバーも同様である。
実際に、補助画像Pは、画像空間またはk空間において従来の方法（例えば、Tihonov正則化を用いて）で非正則化または二次正則化された、同一のR値で１D-SENSE再構成を用いて得ることができる。

それゆえ、最尤（ML）推定手順は、この補助画像が完全事前分布の実現であることを仮定し、それに直接Θをフィットさせることからなる：

統計学において、この解法は、完全データ最大尤度と呼ばれ、前記欠測データ最尤推定量とは対照的である。この手順は、２つの独立したステップに分解でき、第１のステップは、近似係数ζバー_aに付随するガウシアン事前パラメータ(μ,σ)の設定に関し、第２のステップは、GGL事前パラメータの推定に関する。

対応する詳細係数から

一方、ML推定量(μハット,σハット)は、経験平均および標準偏差によって明示的に与えられる

各分解能ｊおよび配向oに対し、

は、以下のようにζバー_o,jから推定される。

ハイパーパラメータ

が、(26)式において、Re(・)をIm(・)に置き換えることによって、同様に推定される。

この二次元最小化問題は、閉形式解を認めない。それゆえ、MLパラメータは、直接探索法（例えば、Rosenbrock法、[Bertsekas、1995、第１章、p159-165]参照）のように数値最適化手法を用いて計算される。モンテカルロ法またはシュタイン方式に基づく代替法は、また、増大する計算負荷を犠牲にしていると考えられている。

時間的正則化を考慮するとき、２つの追加パラメータを決定しなければならない：時間的正則化パラメータκと形状パラメータｐである（21および22を参照）。

実際、被験者レベルの分析において、１つの被験者について記録されたデータのみが処理される。このため、正則化パラメータκおよび形状パラメータｐを手動で決定することが、例えこれが、提案された方法を完全に自動にせず、いくらかユーザに依存したままであるとしても、十分な量のデータについて適切な場合がある。それゆえ、グループレベルの分析を考慮すると、処理データはずっと大きくなり、多くの被験者を考慮すると、一度に全ての被験者について手動で時間的正則化パラメータを設定することは、被験者間のばらつきのため、必ず準最適になる。このようなばらつきは、解剖学的または機能的な観点のいずれにおいても重要なことがある。別々の被験者ごとに、このパラメータを手動で設定することは可能だろうが、煩わしく、提案された方法を（一般に１５の被験者が関係する）グループレベルの分析において関係する被験者の数について一層ユーザに依存させるだろう。

グループレベルの分析においてさえ、別の被験者ごとに推定を実施する。さらに、３Ｄの各ボクセルに対し、パラメータ値をなるべく推定し、このことは、各ボクセルが他のボクセルとは独立に処理されることを意味する。それにも関わらず、計算効率のため、全体積について同一の形状パラメータｐを考慮する。別の言い方をすれば、一意的な時間的正則化パラメータκ、それぞれのパラメータκ_rは、各ボクセルに対して決定される。

前記パラメータの最尤推定を再度実施する。ρ¹(r)…ρ^Nr(r)を、４Ｄデータセットを形成する３Ｄ画像として、V_rをのように定義する：

その後、ＭＬSENSEにおけるパラメータρおよびκ_rの共同推定は、以下の基準を最小化することによって実施される。

好ましくは、ｐ≧１の制約下である。

明らかに、解剖学的および機能的な観点と同様に、所定のボクセルがその隣に似た振る舞いをする可能性が高いと想定される。代替案は、脳の機能領域ごとにパラメータ値を推定することからなる（すなわち、類似の機能的役割を有する隣接するボクセルのグループ）。これは、異なる機能領域を特定するため、脳の事前分類手順に基づく。このような拡張は、より多くのサンプルが与えられた推定手順で利用可能となるため、推定過程をより強固なものにする。実際、同一の機能的領域に属する全ボクセルに関する時間的信号は、同一のベクトル内で集められる。機能的領域ιの大きさをＳιによって示すなら、最尤推定は、Ｓι倍のサンプルに関する信号に基づいてなされる。それゆえ、結果として生じる推定は、ｆＭＲＩセッションにおいて発生することがある、より強固な非定型の観測（異常値）である。

[Chaariら2008]および[Chaariら2009]の画像再構成法は、ウェーブレット基底上の観測体積の分解に基づく。この発明の他の面を構成する、さらなる改善点は、フレーム上で分解が実施される場合に前記方法を一般化することにある。

フレームは、ウェーブレット、カーブレット（curvelet）、バンドレット（bandlet）他を一般化する。フレームは、信号（または単に多次元信号である、画像）が分解できる関数族である。一般に、結果として得られる信号表現は冗長である；冗長でないとき、そのフレームは基底と呼ばれる。フレームは、一般に調査中の画像のより幾何学的な詳細の収集を促進するために用いられる。このため、再構成アーチファクトを追跡するために、ｐＭＲＩ再構成においてフレームを用いることが役立つ。

さらに正確には、それぞれ＜・｜・＞および||・||として示される通常のスカラー積およびノルムユークリッド空間Ｒ^L空間の要素として長さＬの実数値デジタル信号を考慮しよう。

Ｋは、Ｌよりも大きいか、または等しい整数とする。以下のようになるように、]0，+ ∞[中の定数μが存在するとき、有限次元空間Ｒ^K空間のベクトル群(e_k)_1≦k≦Kはフレームである。

不等式（28）が等式になる場合、(e_k)_1≦k≦Kタイトフレームと呼ばれる。有界線形フレーム分析作用素Ｆおよび随伴統合フレーム作用素Ｆ^*は、次のように定義される。

Ｆ^*が全射であるのに対し、Ｆが単射であることに留意する。Ｆ¹＝Ｆ^*であるとき、(e_k)_1≦k≦Kは正規直交基底である。冗長フレームの簡単な例は、正規直交基底の結合である。この場合、フレームはμ＝Ｍでタイト（tight）であり、それゆえＦ^*Ｆ＝ＭＩである。ここで、Ｉは恒等作用素である。

以下において、より複雑さを増した２つのフレーム推定問題に取り組む、１つめは、ノイズ除去の定式化に対応する問題であり、２つめは、平衡ＭＲＩ再構成において感度行列Ｓのような分解作用素に関する拡張に対応する問題である。最初にノイズ除去問題に取り組み、そのｐＭＲＩへの応用は、以下ではｙとして示されるノイズの多い基準画像（例えば、ＳＥＮＳＥ再構成のような）がフレーム表示のハイパーパラメータを推定するのに役立つ。第２の問題は、アルゴリズム４を示す際に後述する。ｐＭＲＩの状況で縮小されたＦＯＶ画像から直接フレーム分解のハイパーパラメータを推定することに対応する。

観測された信号y=F^*x+nは、以下のようにx∈Ｒ^K空間に関する係数フレーム表現（ＦＲ）によって記載できる。

ここで、ｎは、観測された信号とＦＲ Ｆ^*ｘとの間の誤差である。この誤差は、それが閉凸集合に属することを課することによってモデル化される。

ここで、δ∈[0, ∞[ は、誤差の境界であり、ＮはＲ^K空間上の任意のノルムである。

信号／画像回復問題において、ｎは測定データを破損する付加雑音である。以下の展開は、一様なノイズによってモデル化された境界観測誤差の場合に焦点を当てる。確率論的アプローチを採用することによって、ｙおよびｘがランダムベクトルＹおよびＸを実現したものと想定される。このような状況において、目的は、パラメトリック確率論的モデルを考慮し、関連するハイパーパラメータを推定することによって、Ｘ｜Ｙの確率分布を特徴付けることである。推定問題は、非冗長基底、例えばウェーブレット基底上での分解の場合においてずっと困難である。なぜなら、Ｆ^*は全単射ではないからである。

ベイズ的な枠組みにおいて、フレーム係数のための事前分布を定義する必要がある。例えば、この事前は、表現の希薄を促進するために選択することがある。以下のf(x, θ)において、未知のハイパーパラメータベクトルθおよびｆ（θ）に依存し、ハイパーパラメータベクトルθ用の事前ｐｄｆであるフレーム係数の確率密度関数（ｐｄｆ）を示す。制約C_δに従って、ｎは、球上に不均一に分布しているものと想定される。

f(y|x) は、凸の閉じた球上の均一なｐｄｆである。

ハイパーパラメータベクトルθはまたランダムベクトルΘの実現も考慮しており、(X, Θ)の条件付きｐｄｆは、以下のように書かれる：

簡単のため、上記のように例え、ＧＧＬ分布がｐＭＲＩへの適用により適しているとしても、フレーム係数が周辺ＧＧ（一般化ガウシアン）分布について事前独立であると仮定する。これは、以下のフレーム係数を事前に生じる：

ここで、α_k>0, β_k>0 (k∈[1, …, K])は、フレーム係数ベクトルの要素であるX_kに関連する尺度および形状パラメータであり、Γ（・）はガンマ関数である。

を導入することにより、フレーム事前分布は以下のように書くことができる。

フレーム係数の分布は一般に互いの係数ごとに異なる。しかしながら、いくつかのフレーム係数は、非常によく似た分布を有することができ、同一のハイパーパラメータによって定義できる。結果として、フレーム係数をＧの異なるグループに分割することが提案される。ｇ番目のグループは以下のように示される一意のハイパーパラメータベクトルによってパラメータ化される。

この場合において、フレーム事前分布は、以下のように表すことができる。

ここで、和は、n_g要素およびθ_g=(θ₁…θ_G)を含むg番目のグループの要素の添数集合Ｓ_gの範囲にわたる。例えば、各グループは所定のウェーブレットサブバンドに対応できる。所定の分解能の全てのフレーム係数が同一のグループに属することを考慮することによって、マルチスケール表示を用いるとき、より粗い分類が多くある場合がある。
フレーム分解のための階層ベイズモデルは、以下の不適切なハイパー事前分布によって完成される：

ここで、ξ∈Aの場合、1_A(ξ)=0、それ以外の場合、1_A(ξ)=1とすることにより、1_A(ξ) は、A⊂Ｒ空間上で定義される関数である。
この種の事前分布を用いることの動機付けは以下の通りである。
・区間［０，３］は実際の応用において遭遇するβ_gの全ての起こりうる値の範囲であり、そのパラメータについて追加の情報はない。
・パラメータγ_gのための事前分布は、このパラメータについての知識の欠如を反映するジェフリー分布である。
結果として得られる事後分布はそれゆえ次式によって与えられる：

事後分布（35）に関連するベイズ推定量［例えば、最大事後（ＭＡＰ）または最大平均二乗誤差（ＭＭＳＥ）推定量も事後平均推定量として知られている］は、単純な閉形式表現を有しない。

以下において、ハイブリッドマルコフ連鎖モンテカルロ（ＭＣＭＣ）アルゴリズムに依存する、ＭＭＳＥ推定量を計算するために確率論的な手順が提案されている。この考えは、事後分布（35）に従って分布するサンプルを生成するために適したアルゴリズムを用いることである。収束後、生成されたサンプルが、それぞれ未知のモデルパラメータおよびハイパーパラメータベクトルＸおよびθのＭＭＳＥ推定量を計算するために用いられる。ｐＭＲＩの状況において、ｘはリファレンス画像の得られたフレーム表現である。

考慮されたフレームが正規直交Ｍ基底およびＮ（・）がユークリッドノルムであるとき、対象の分布と関連する条件分布に従って分布したサンプルを繰り返し生成するよく知られたギブスサンプリング法（ＧＳ）を用いることができる［Ｇｅｍａｎ、１９８４］。さらに正確には、基本的ＧＳは完全事後f(x,θ|y)の実現をシミュレートするf(x|θ,y)およびf(θ|x,y)に従って分布したサンプルを繰り返し生成する。

直接計算によって、以下の条件分布が得られる。

この条件分布は、C_δ上で切断されたＧＧ分布の積である。実際、この切断された分布によるサンプリングは、必ずしも実施するのが容易ではない。なぜなら、随伴フレーム作用素Ｆ_は通常大きな次元であるからである。しかしながら、２つの代替サンプリング戦略が以下のように詳細に説明される。

単純サンプリングは、独立したＧＧ分布

に従うサンプリングによって進められる。その後、N(y-F^*x)≦δの場合のみ、提案された被験者Ｘを採択する。この方法は、任意のフレーム分解および任意のノルムに対して用いることができる。しかしながら、それは全く非効率である。なぜなら、特にδが小さい値をとるとき、採択比率が非常に低いためである。

ギブスサンプリング法は、考慮されたフレームがＭ個の正規直交基底であり、Ｎ(・)がユークリッドノルムであるとき、（36）式の条件分布から、より効率的なサンプルに設計される。この場合において、解析フレーム作用素および対応する随伴行列は、それぞれF=[F₁…F_M]^T および F^*=[F^* ₁…F^* _M]^Tとかくことができる。ここで、∀m∈[1,…,M], F_mは、F^* _mF_m=F_mF^* _m=Idのようなｍ番目の正規直交基底上への分解作用素である。以下において、Ｋ＝ＭＬの各x∈Ｒ^K空間は、x=[x₁ ^T,…,x_M ^T]^Tのように分解される。ここで、x_m∈Ｒ^L ∀m∈[1,…,M]である。

フレーム係数の生成用のＧＳは、制限N(y-F^*x)≦δのもとでの条件分布f(x_n|x_-n,y,θ)によるベクトルをひく。ここで、x_-nは、ｎ番目のベクトルx_nを取り除くことによって、ｘから形成されたＲ^K-L空間の次元が縮小されたサイズのベクトルである。N(・)がユークリッドノルムの場合、∀n[1,...,M]は:

である。
各x_nをサンプルするため、提案された分布が以下の式によって定義される球B_Cn,δ上で支持されるＭＨステップを用いることが提案されている。

上で定義されるｐｄｆ q_δ のランダム生成がこれから議論されるだろう。

最初に、Ｒ^L空間の単位l_p（p ∈]0,+∞]）の球内のベクトルをどのようにサンプリングすべきかが考慮される。特別な場合ｐ＝＋∞において、区間[-1,1]上の分布に従って各空間座標に沿って独立にサンプリングすることによって、これを容易に実施できる。ここで考慮されたパラメータ“ｐ”は、時間的な正則化で用いられる形状パラメータと混同すべきでない。21式および22式を参照せよ。
このように、このセクションは、ｐの有限値と関連する一層困難な問題に焦点を当てる。以下において、||・||_pはｃ_pノルムを示す。
A=[A₁,…,A_L']^Tを以下のGG(p^1/P,p)pdfを有するi.i.d.要素の確率ベクトルとする。

U=[U₁, … , U_L']^T = A/||A||_pとする。その後、乱数ベクトルＵは、Ｒ^L'空間のｌ_p単位球の表面およびU₁,…,U_L'-1の結合ｐｄｆの一様に分布していることを示すことができる（[Guptaら1997]）。

Ｒ^L'空間のｃ_p単位球の一様な分布は、u(L',p)によって示される。
それゆえ、U=[U₁, … ,U_L']^T 〜 u(L',p)。
各[1,…,L'-1]に対して、V=[U₁, … ,U_L]^T のpdfは次式によって与えられることを示すことができる([Guptaら1997])：

特に、ｐ∈Ｎ^*空間およびL'=L+pなら、(40)式はＲ^L'空間のｃ_p単位球の一様な分布を提供する。半径l_pの球（ball）上の任意の分布ｑ_nからサンプリングして、

が、Ｖを推定することによって直接推定できる。

このようなｐｄｆの閉形式表現を有することは、ＭＨ生成の採択率を計算することができるため重要である。前の反復（i-1）で得られたx_n ^(i-1)の値を考慮に入れるために、x_n ^(i-1)を含む半径η∈]0,δ[の限定された球上に支持される提案された分布を選択することが好ましいことがある。ランダムウォークＭＨに似たこの戦略は、条件分布ｆ(x|θ,y)の大きな値に関連する領域のより良い調査をもたらす。より正確には、次式において定義される提案分布を選択することが提案される。

そして、Ｐは以下のように定義される球B_0,δ-η上への射影である。

このような球の中心の選択は、以下の関係を保証する。

さらに、以下のいかなる点も、

内での連続的な引用後に到達できる。

以下をよく調査を確保するために、半径ηを調整しなければならない。

実際、採択率を改善するために（δと比べて）ηの十分に小さい値を決めることが興味深いこともある。
f(θ|x, y)によってハイパーパラメータベクトルθをサンプリングする代わりに、f(γ_g|β_g,x, y) および f(β_g|γ_g,x, y)によって繰り返しサンプリングすることが提案されている。直接的な計算によって、以下の結果が得られる：

結果として、f(γ_g|β_g,x, y)は、サンプリングが容易な逆ガンマ分布のｐｄｆである。

反対に、切断されたｐｄｆ f(β_g|γ_g, x, y)によってサンプリングすることはより困難である。これは、提案されたq(β_g|β_g ^(i-1))が標準偏差δβ_g= 0.05を有し、区間[0,3]で切断されたガウス分布である、ＭＨ生成を用いることによって実現できる。この分布モードは、前の反復(i-1)でのパラメータβ_g ^(i-1)の値である。結果として生じる方法は、アルゴリズム４において要約される合成ＧＳである。
もし、同一のベイズ構成内での再構成およびハイパーパラメータ推定手順を１つにまとめたいのなら、30式の観測モデルをy=SF^*x+nに拡張できる。ここで、Ｓは、(3)式で定義される感度行列である。その結果、推定が再度実施されねばならない。

アルゴリズム４
１：θ⁽⁰⁾＝（θ_g ⁽⁰⁾）_1≦g≦G＝（γ_g ⁽⁰⁾、β_g ⁽⁰⁾）_1≦g≦Gおよびｘ⁽⁰⁾∈Ｃ_δで開始してi=1と設定
２：repeat
３：ｘをサンプリングする：
４：for 1 to M do
５：

およびｘハット_n ^(i-1)＝Ｐ（ｘ_n ^(i-1)−ｃ_n ⁽ⁱ⁾）＋ｃ_n ⁽ⁱ⁾を計算
６：ｘ_n ⁽ⁱ⁾を以下のようにシミュレートする：
・ｘチルダ_n ⁽ⁱ⁾〜ｑ_η（ｘ_n−ｘハット_n ^(i-1)）
ここで、ｑ_ηはＢ_0,η上で定義される。
・比率

を計算し、確率ｍｉｎ｛１，ｒ（ｘチルダ_n ⁽ⁱ⁾，ｘ_n ^(i-1)）｝で提案された被験者を採択する。
７：end for
８：θをサンプリングする：
９：for g=1 to G do
１０：

を生成
１１：以下のようにβ_g ⁽ⁱ⁾をシミュレートする。
・βチルダ_g ⁽ⁱ⁾〜ｑ（β_g｜β_g ^(i-1)）
・比率

を計算し、確率ｍｉｎ｛１，ｒ（βチルダ_g ⁽ⁱ⁾，β_g ^(i-1)）｝で提案された被験者を採択する。
１２： end for
１３：ｉ←ｉ＋１
１４：until 収束

このアルゴリズムは、直感で理解でき、実施するのが容易であるが、考慮されたフレームがＭ個の正規直交基底の結合であるという限定的な仮定の下に導かれたということを指摘しなければならない。これらの仮定が後述の他のアルゴリズムに当てはまらないとき、フレームの代数的な性質を利用することによってフレーム係数および関連するハイパーパラメータのサンプリングを可能にする。実際、f(x|θ,y)によるサンプリングの直接的生成は、一般に不可能であるため、ギブス生成をＭＨ生成に代わる代替方法が提案されている。このＭＨ生成は、提案分布による被験者ｘを全体的にサンプリングすることを目的とする。この被験者は、標準的なＭＨ採択率で採択または棄却される。ＭＨ生成の効率は、ｘについての提案された分布の選択に強く依存する。

ｘ（ｉ）をアルゴリズムのｉ番目に採択されたサンプルであり、q(x|x^(i-1))を反復ｉで対象を生成するために用いられる提案とする。q(x|x^(i-1))を選択するための主な困難は、q(x^(i-1)|x)/q(X|X^(I-1))の扱いやすい表現を生じる一方で、x∈C_δを保証しなければならないという事実から生じる。

このため、フレーム表現の代数的な性質を利用することが提案されている。より正確には、いかなるフレーム係数ベクトルもx=x_H+x_H⊥として分解できる。ここで、x_Hおよびx_H⊥はそれぞれ、H = Ran(F)およびH^⊥, = [Ran(F)]^⊥=Null(F^*)において、それらの値をとる確率ベクトルが具体化したものである。Fの範囲は、Ran(F)=[x∈Ｒ^K空間|∃y∈Ｒ^L空間, Fy=x]であり、Ｆ^*の零空間は、Null(F^*)=[x∈Ｒ^K空間|F^*x=0]であることを想起する。

ここで用いられる提案された分布は、サンプルx_H∈H および x _H⊥∈H^⊥を生成することを可能にする。より正確には、提案されたｐｄｆの以下の分離形式が考慮される：

ここで、x_H ^(i-1)∈H, x_H⊥ ^(i-1) ∈H^⊥ および x^(i-1)= x_H ^(i-1) + x_H⊥ ^(i-1)である。言い換えれば、x_H および x_H⊥の独立サンプリングが実施される。

分解x =x_H+x_H⊥を考慮すると、Ｃ_δにおいてｘをサンプリングすることは、λ∈Cバー_δのサンプリングに等しい。ここで、Cバー_δ=[λ∈Ｒ^L空間|N(y-F^*Fλ)≦δ]である。実際、x_H=Fλwであり、ここで、λ∈Ｒ^L空間であり、そしてx_H⊥∈Null(F^*)からF^*x=F^*Fλである。Ｃバー_δにおけるλのサンプリングは、例えば、球B_y,δ上の分布からｕを生成することによって、そしてλ=(F^*F)^-1uととることによって容易に実施できる。

反復ｉでのｘＨのサンプリングをより効率的にするために、前の反復x_H ^(i-1)=Fλ^(i-1)=F(F^*F)^-1u^(i-1)のサンプリング値を考慮することが興味深い場合がある。

ランダムウォーク生成技法と同様に、

において、ｕが生成される。ここで、η∈]0, δ[ および uハット^(i-1)=P(u^(i-1)-y)+yである。これは、x_H=F(F^*F)^-1u∈C_δ および N(u-u^(i-1))≦2ηのようなベクトルｕを描くことを可能にする。
その後、N(・)は、p∈[1,+∞]を有するc_p ノルムであると仮定して、ｕの生成が(40)式により実施できる。

一度、x_H=Fλ ∈ H∩C_δ (ｘはC_δ内にあることを保証する)がシミュレートされると、H_⊥の要素としてサンプリングしなければならない。y=F^*x+n=F^*x_H+nであるため、x_H⊥についてｙに情報がない。結果として、ガウス分布N(x^(i-1),σ_x ²I)に従ってｚを描き、H^⊥上にｚを射影すること、すなわち、x_H=П_H⊥zによって、x_H bをサンプリングすることが提案されている。

それゆえ、提案ｐｄｆの表現は、以下の式であることを示すことができる。

この表現は、Ｋ＝Ｌ（x _H⊥=0を生じる）であるとき、縮退の場合において有効のままである。最後に、ｑ_ηが球B_{0, η}上の一様な分布のとき、前記比率はＭＨ採択比率を簡単にする１に帰着することに注目することが重要である。最後のアルゴリズムは、アルゴリズム５に要約される。アルゴリズム４のハイブリッドＧＳについて、ハイパーパラメータベクトルのサンプリングを実施することに注目すべきである。

アルゴリズム５
１：θ⁽⁰⁾＝（θ_g ⁽⁰⁾）_1≦g≦G＝（γ_g ⁽⁰⁾、β_g ⁽⁰⁾）_1≦g≦Gおよびｕ⁽⁰⁾∈Ｂ_y,δで開始
ｘ⁽⁰⁾＝Ｆ（Ｆ^*Ｆ）^-1ｕ⁽⁰⁾とおき、i=1と設定
２：repeat
３： ｘをサンプリングする：
・ｕハット^(i-1)＝Ｐ（ｕ^(i-1)−ｙ）＋ｙを計算
・ｕチルダ⁽ⁱ⁾〜ｑ_η（ｕ−ｕハット^(i-1)）を生成
ここでｑ_ηは、Ｂ_0,η上で定義される。
・ｘチルダ_H ⁽ⁱ⁾＝Ｆ（Ｆ^*Ｆ）^-1ｕチルダ⁽ⁱ⁾を計算
・ｚ⁽ⁱ⁾〜Ν（ｘ^(i-1)，σ_x ²Ｉ）を生成
・ｘチルダ_H⊥ ⁽ⁱ⁾＝Π_H⊥ｚ ⁽ⁱ⁾および
ｘチルダ⁽ⁱ⁾＝ｘチルダ_H ⁽ⁱ⁾＋ｘチルダ_H⊥ ⁽ⁱ⁾を計算
・比率

を計算し、提案された対象ｕチルダ⁽ⁱ⁾およびｘチルダ⁽ⁱ⁾を確率ｍｉｎ｛１，ｒ（ｘチルダ⁽ⁱ⁾，ｘ^(i-1)）｝で採択する。
４： θをサンプリングする：
５： for ｇ＝１ to Ｇ do
６：

を生成する
７： 以下のようにβ_g ⁽ⁱ⁾をシミュレートする
・βチルダ_g ⁽ⁱ⁾〜ｑ（β_g｜β_g ^(i-1)）を生成する
・比率

を計算し、確率ｍｉｎ｛１，ｒ（βチルダ_g ⁽ⁱ⁾，β_g ^(i-1)）｝で採択する。
８： end for
９： ｉ←ｉ＋１
１０：until 収束

再構成すべき画像の全変分（ＴＶ）に依存する事前分布における追加項を含むように、（31-35）式の階層ベイズモデルを拡張できる。冗長フレーム表現のように、ＴＶ事前分布を用いることがハイパーパラメータ推定問題をも生じる。
指数関数型を想定して、新しい事前分布は以下のように表現できる：

ここで、θ=(θ₁, … , θ_G, κ)は、κ＞０の場合の新しいハイパーパラメータベクトルであり、||・||_TV はＴＶ半ノルムであり、Z(θ)は正則化定数である。フレーム分解のための新しい階層ベイズモデルは、以下の不適正ハイパー事前分布によって完結する：

ここで、κ_maxは、（好ましい実施例においては、κ_max＝１０に）固定すべき正の実数である。

Ｚ（θ）は、γ_gに関して一様に有界であり、それゆえ(47)式のハイパー事前分布ｆ（θ）は、γ_g→+∞のとき安定漸近特性を有することを示すことができる。
その結果生じる新しい事後分布は、それゆえ以下の式によって与えられる：

(48)式の事後分布に関連したベイズ推定量はなお、単純な閉形式の表現を有しない。この理由のため、上述したものと同一のサンプリング戦略が適用され、フレーム係数がアルゴリズム５のようにサンプリングされる。しかしながら、ハイパーパラメータベクトルについては、直接計算によって、ハイパーパラメータγ_g β_g およびκの事後分布がそれぞれ次式のように表現される：

結果として、f(γ_g|β_g,κ,x,y)は、γの逆数のｐｄｆである。

γ_gのサンプリングは、それゆえアルゴリズム４のように正確に実施される。
逆に、f(γ_g|β_g,κ, x, y)およびf(γ₁, … ,γ_G,β₁, …, β_G, x, y)によるサンプリングはより困難である。このような課題は、提案された分布q(β_g|β_g ^(i-1))およびq(κ|κ^(i-1))がそれぞれ、区間［０，３］および[0, κ_max]、σβ_g=0.05およびσκ =0.01の標準偏差を有するガウス分布で切断された２つのＭＨ生成を用いることによって実現できる。これらの標準偏差値は、実際の観測結果に基づき決定される。(49)式の事後分布によってサンプリングするために結果として生じる方法は、アルゴリズム６に要約される。

アルゴリズム６
１：θ⁽⁰⁾＝（（θ_g ⁽⁰⁾）_1≦g≦G、κ⁽⁰⁾）＝（（γ_g ⁽⁰⁾、β_g ⁽⁰⁾）_1≦g≦G、κ⁽⁰⁾）およびｕ⁽⁰⁾∈Ｂ_y,δで開始
ｘ⁽⁰⁾＝Ｆ（Ｆ^*Ｆ）^-1ｕ⁽⁰⁾とおき、i=1と設定
２：repeat
３： ｘをサンプリングする：
・ｕハット^(i-1)＝Ｐ（ｕ^(i-1)−ｙ）＋ｙを計算
・ｕチルダ⁽ⁱ⁾〜ｑ_η（ｕ−ｕハット^(i-1)）を生成
ここでｑ_ηは、Ｂ_0,η上で定義される。
・ｘチルダ_H ⁽ⁱ⁾＝Ｆ（Ｆ^*Ｆ）^-1ｕチルダ⁽ⁱ⁾を計算
・ｚ⁽ⁱ⁾〜Ν（ｘ^(i-1)，σ_x ²Ｉ）を生成
・ｘチルダ_H⊥ ⁽ⁱ⁾＝Π_H⊥ｚ⁽ⁱ⁾および
ｘチルダ⁽ⁱ⁾＝ｘチルダ_H ⁽ⁱ⁾＋ｘチルダ_H⊥ ⁽ⁱ⁾を計算
・比率

を計算し、提案された対象ｕチルダ⁽ⁱ⁾およびｘチルダ⁽ⁱ⁾を確率ｍｉｎ｛１，ｒ（ｘチルダ⁽ⁱ⁾，ｘ^(i-1)）｝で採択する。
４： θをサンプリングする：
５： for ｇ＝１ to Ｇ do
６：

を生成する
７： 以下のようにβ_g ⁽ⁱ⁾をシミュレートする
・βチルダ_g ⁽ⁱ⁾〜ｑ（β_g｜β_g ^(i-1)）を生成する
・比率

を計算し、確率ｍｉｎ｛１，ｒ（βチルダ_g ⁽ⁱ⁾，β_g ^(i-1)）｝で採択する。
８： end for
９： 以下のようにκ⁽ⁱ⁾をシミュレートする
・κチルダ⁽ⁱ⁾〜ｑ（κ｜κ^(i-1)）を生成する
・比率

を計算し、提案された対象を確率ｍｉｎ｛１，ｒ（κチルダ⁽ⁱ⁾，κ^(i-1)）｝で採択する。
１０：ｉ←ｉ＋１
１１：until 収束

この発明の方法の様々な実施形態の技術的結果は、図２〜２０を参照することによってこれから議論される。

図２〜１０および１３で報告された解剖学的な結果は、以下の条件で得られたものである。０．９３×０．９３×８ｍｍ^３の空間分解能で解剖する２５６×２５６×１４のグラジェントエコー（ＧＥ）を備えた実データセットで実験が行われた。ＧＥ解剖学的画像が、TE/TR=10/500msおよびBW=31.25kHzで収集された。これらの画像は、８チャネルヘッドコイルを有するSigna １．５テスラＧＥヘルスケアスキャナで短縮ファクタＲ＝２およびＲ＝４を用いて収集されたということも留意すべきである。興味深いことに、解剖学的データの走査時間は、非パラレル撮像法において５〜ｍｎ続いたものの、パラレル撮像法においてＲ＝２およびＲ＝４でそれぞれ収集時間が３ｍｎ１０ｓおよび２ｍｎ２０ｓに減少した。

図２および３は、ＳＥＮＳＥを用いてＴｉｋｈｏｎｏｖ正則化で得られた人間の脳の９枚の解剖学的スライスを示す。Ｔｉｋｈｏｎｏｖパラメータは手動で設定されている。減少係数は、図２でＲ＝２、そして図３でＲ＝４である。正則化に関わらず、いくつかのエイリアシングアーチファクトが、特に、図３上において曲線形状において見られる。画像の過剰な平滑化も見られる。

同一の実験条件でＴＶ正則化（左側のパネル：Ｒ＝２；右側のパネル：Ｒ＝４-単一のスライスが示されている）を用いて、手動で調整された正則化パラメータκで図４が得られた。かなり強いエイリアシングアーチファクトおよび“階段”欠陥が見られる；これらの欠陥は、画像の情報量を犠牲にしてκ値を増加させることによって低減できる。

上述のように、同じ実験条件で、自動較正された２Ｄのウェーブレット変換基底に基づく正則化方式を用いて、図５および６が得られた；この手法は、２Ｄ−ＵＷＲ−ＳＥＮＳＥと呼ばれるが、“ＵＷＲ”は“制約なしのウェーブレット正則化”を表す。また、図５でＲ＝２、そして図６でＲ＝４である。

より正確には、長さ８のフィルタに関連する２進（Ｍ＝２）Ｓｙｍｍｌｅｔ正規直交ウェーブレット基底がｊ_max＝３の分解能以上で用いられた。ウェーブレット係数について、（10）式の事前分布が採用されている。

関連するハイパーパラメータ（一組のハイパーパラメータが、各サブバンド、すなわち各分解能および方向での各近似／詳細係数の実／虚部についてフィットされる）が、上述のベイズ的アプローチを用いて推定された。その後、完全ＦＯＶ画像の再構成がウェーブレット基底に基づく正則化を用いて実施された。

Ｔｉｋｈｏｎｏｖ正則化について図２および図３で観測された平滑作用は、図５および６のＷＴ正則化画像においてもはや存在しない。ここで、ＴＶ正則化に関して観測される階段効果を導入することなしに、かなり正確な再構成が脳マスク内で実施される（図４参照）。

アーチファクト領域をより良く正則化するため、上述の方法において追加の制約を導入することによって画質のさらなる改善が得られる。次のアルゴリズムは、ＣＷＲ−ＳＥＮＳＥと呼ばれており、“ＣＷＲ”は“制約付きウェーブレット正則化”の略称である。

この方法は、アーチファクト領域の形状および／または位置に関わらず、アーチファクト領域内の画像強度値の局所的な下限および上限を課すことを意味する。これらの境界は、非空凸閉包集合を規定する：

ここで、位置r∈2 [1, … ,Y/R] × [1, … ,X]の範囲の値で導入される制約は、以下によってモデル化される：

最適性基準は次のようになる：

ここで、C^*=[ζ∈Ｃ空間^K|T^*ζ∈C] およびiC^* は、次式によって定義される閉包集合Ｃ^*の指標関数である：

追加の凸制約が適用されるアーチファクト領域を見つけるため、（勾配画像において非常に低い／高い推移の）形態学的勾配[Serra, 1982]を用いることができる。凸集合Ｃ_rを規定する上限および下限は、ｒに依存するため空間的に変化する。それゆえ、これらは非常に低く高い強度を切り捨てるため、基本ＳＥＮＳＥ再構成画像（補助またはリファレンス画像）に適用される形態学的開閉操作を用いて計算できる。

図７および８は、ＣＷＲ−ＳＥＮＳＥ法の結果を、それぞれＲ＝２およびＲ＝４について例示する。追加の凸制約を用いて異方的な平滑化によって、図５および６に残存しているアーチファクトが現在取り除かれていることを見ることができる。

定量的な観点から、基本ＳＥＮＳＥおよびＴｉｋｈｏｎｏｖ再構成と比べて、ＵＷＲ／ＣＷＲ−ＳＥＮＳＥによって著しい改善が達成された。以下の表１は、図３，６および８（Ｒ＝４）において示された解剖学的な脳の体積の例示されたスライスについての基本−ＳＥＮＳＥ、Ｔｉｋｈｏｎｏｖ正則化および提案されたＵＷＲ／ＣＷＲ−ＳＥＮＳＥに対応する信号対雑音比（ＳＮ比）のｄＢ値を報告する。提案された制約付きの正則化戦略から得られたゲインは、平均して１．０８ｄＢになり、より良い視覚品質である。

ウェーブレット基底の選択の影響も研究されている。特に、４つの異なる基底が考慮されている：Ｍ＝４バンドを有する２進 Symmlet 8, ２進 Daubechies 8, ２進 Haar および Meyerである[Chauxら2006b]。最初の３つの基底は（Ｈａａｒ基底でのみ生じる何らかの遮断効果は別として）かなり似た結果を与え、２進 Symmlet 8はわずかに高いＳＮ比を生じる。それどころか、Ｍｅｙｅｒ４バンドウェーブレット基底は、著しく低いＳＮ比を生じる。

図９は、Ｊ_max＝３以上の分解能を有する長さ８フィルタに関連する２進（Ｍ＝２）Ｓｙｍｍｌｅｔ正規直交ウェーブレット基底を減少係数Ｒ＝４で用いた結合ウェーブレット全変分正則化（アルゴリズム２）によって得られた結果を示す。図１０は、Ｓｙｍｍｌｅｔ８およびＳｙｍｍｌｅｔ４のフィルタを有する２つの正規直交基底の結合（Ｕ２ＯＢ）によって構成された冗長ウェーブレットフレームを用いた結合ウェーブレット全変分正則化によって得られた結果を示す。ウェーブレット事前分布のハイパーパラメータおよびＴＶ正則化パラメータが上記手法を用いて推定された。これらの図は、再構成された画像がＵＷＲ−ＳＥＮＳＥアルゴリズムを用いて再構成されたものよりもより良い規則性を提示することを示す（図６）。しかしながら、視覚的な観点から冗長ＷＴを用いることは必ずしもより良い再構成の品質を生じるものではない。定量的な観点から、ＳＮ比のｄＢ値が以下の表２に示される。表１のＳＮ比値と比較すると、結合正則化構造においてＷＴとＴＶを結合することの有用性が裏付けられる。この表は、Ｓｙｍｍｌｅｔウェーブレットを用いるとき、ＵＷＲ−ＳＥＮＳＥアルゴリズムと比べて、０．０２ＤＢのわずかな改善が得られることを示す。しかしながら、Ｕ２ＯＢ冗長ウェーブレットフレームを用いるとき、より顕著な改善（０．２４ＤＢ）が達成される。たとえ正則化ウェーブレット基底および冗長ウェーブレットフレームを用いて得られる再構成性能が等しいように思われたとしても、冗長フレームを用いることが有益であることを、ＳＮ比値が示すことが、その後判明している。

図１３は、ＴＶ正則化（左列）、２つの正規直交ウェーブレット基底の結合によって構成されたウェーブレットフレームを用いる３Ｄ−ＵＷＲ−ＳＥＮＳＥ（中列）およびハイブリッドフレームＴＶ正則化法（右列）を用いた再構成画像を示す。上段は全体画像を示し、下段は拡大詳細図を示す。その結合が再構成用に用いられるウェーブレットフレームを構成する２つの正規直交基底（中列および右列）は、長さ４のＤａｕｂｅｃｈｉｅｓ規定および長さ８の変換Ｏａｕｂｅｃｈｉｅｓ規定である。ウェーブレット係数のＧ＝２０グループが考慮されることを意味する、３つの分解能が用いられている。フレーム表現のため、上述のハイブリッドＭＣＭＣアルゴリズムを用いてハイパーパラメータが選択されている。

図１１および図１２で報告された解剖学的結果が以下の条件で得られた。３Ｄ Ｔ１強調ＭＰ−ＲＡＧＥパルス系列と３２個の受信チャネルからなるマトリックスアレイヘッドコイルを用いて３Ｔ Ｔｉｍトリオジーメンススキャナで解剖学的ＭＲＩスキャンが実施された。スキャンパラメータは以下のように選択された：スライス方向＝矢状断（右->左）、スライスの厚み＝１．１ｍｍ；ＴＥ＝２．９８ｍｓ、ＴＲ＝２３００ｍｓ；Ｔｌ＝９００ｍｓ、フリップ角：９°、ＢＷ＝６１ｋＨｚ、単発。視界は、１×１×１．１ｍｍ³の異方的分解能に対応する２５６×２４０×１６０のマトリックスサイズを有する２５６×２４０×１７６ｍｍ³であった。走査時間は、従来のＭＲＩを用いて、すなわち短縮なし（ｐＭＲＩでも部分フーリエのいずれでもない）で、ＴＡ_ｃｏｎｖ＝９分１４秒だった。６／８部分フーリエスキームの利用により、走査時間が７分４６秒まで減少させることが可能になった。パラレルＭＲＩのデータが短縮ファクタＲ＝２またはＲ＝４を用いて部分フーリエなしに収集された。これらの試験のための各スキャン時間は、５分０３秒および２分５９秒であり、予想されたようなＴＡ_ｃｏｎｖではなかった。この理由は、ジーメンスのk空間サンプリング戦略にある。製造者は、規定のものより低い実際の短縮ファクタを含む２４本の中心線に対する完全なサンプリングスキームを採用する（２の代わりに１．８３および４の代わりに３．０９）。

図１１〜１２は、マトリックスアレイ３２チャネルヘッドコイルを用いて、３テスラ（ジーメンスＴｉｍトリオ）で収集された単一被験者のＴ１強調ＭＲＩデータからの別のＭＲＩ再構成アルゴリズムを例示する。図１１および１２は、同一の被験者についてそれぞれ軸位断および冠状断像を示す。これらの図の上段は根拠となる事実、すなわち完全なｋ空間の収集から実施される再構成を示す。次に、ｍＳＥＮＳＥ（中段）と３Ｄ−ＵＷＲ−ＳＥＮＳＥ（下段）のアルゴリズムが雑音の多い状況、すなわちＲ＝４について比較される。後者のアルゴリズムについて、Ｊ＝２分解能のレベルおよびＤａｕｂｅｃｈｉｅｓウェーブレット（非冗長）についての結果が得られている。また、ウェーブレット表示のハイパーパラメータが、本明細書において上で詳細に述べた、完全データ最尤手順を用いて推定された。灰白質の境界上の白質のスポットとして現れるＭＳＥＮＳＥ再構成アーチファクトは、いずれの表示の輪郭においても白丸によって示される。これらのアーチファクトは同一位置に現れないため、３Ｄ−ＵＷＲ−ＳＥＮＳＥアルゴリズムは、明らかにｍＳＥＮＳＥ技術よりも優れている。また、軸位断像の中心の楕円形状のアーチファクトはこの発明のウェーブレット基底のアルゴリズムによって強くフィルタされる。冠状断像は、提案された発明を用いて、脳幹および皮質下領域内のかなりの騒音が低減することを裏付ける。

２Ｄウェーブレット基底を用いて得られた結果は、図１１および１２において説明されたものに定性的に近いため、ここでは示さない。しかしながら、３Ｄウェーブレット基底の再構成は、定量的には信号対雑音比（ＳＮ比）の観点から２Ｄウェーブレット基底の再構成よりも優れている。３Ｄ再構成が、２Ｄ技術に対してＳＮ比を１．３ＤＢ改善したことに留意すべきである。

ここまで、解剖学的な静止画撮像法についての結果のみが議論されている。しかしながら、上述のように、例えばエコープラナー撮像法（ＥＰＩ）によって得られた、４Ｄ ｆＭＲＩ画像系列の空間−時間正則化にもこの発明は適用される。対応する結果が、図１４−２０によって例示される。

検証のため、認識ローカライザー[Pinelら2007] プロトコルの間、グラジェントエコーＥＰＩ（ＧＥ−ＥＰＩ）列(TE = 30ms、TR = 2400 ms、スライスの厚み = 3 mm、横断方向、FOV =192 mm²)を用いて、３Ｔジーメンストリオ磁石でｆＭＲＩデータを収集した。この実験は、聴覚、視覚および運動脳機能の地図を作るほか、数の処理および言語理解（リスニングとリーディング）のような、より高い認知的作業も目的としている。それは、Ｎ_r＝１２８スキャンの単一セッションからなる。理論的枠組みは、１０の実験条件（聴覚および視覚文、聴覚および視覚計算、左／右聴覚および視覚クリックで定義された６０個の聴覚、視覚および運動刺激を備えた高速事象に関連するデザインであった。Ｌ＝３２チャネルコイルがパラレル撮像法を可能とするために用いられた。

2×2 mm²の空間面内の分解能で、異なる減少係数（Ｒ＝２またはＲ＝４）を用いて１５の被験者のそれぞれに対してｆＭＲＩデータが収集された。スキャナによって供給された生データファイルに基づき、縮小されたＦＯＶ ＥＰＩ画像が２つの特殊な処理を用いて再構成された。
i）ＧＥ−ＥＰＩのような高速ＭＲＩ列で発生する、読み出し勾配ランプ中に不均一なｋ空間サンプリングを解消するためｋ空間を再構成（regridding）すること。
ii）ＥＰＩ画像のｋ空間収集中に奇数−偶数エコーの不一致を除去するためのナイキストゴースト補正。

その後、４Ｄ−ＵＷＲ−ＳＥＮＳＥおよび既に議論された３Ｄ−ＵＷＲ−ＳＥＮＳＥ（単に“ＵＷＲ−ＳＥＮＳＥ”と略す）アルゴリズムが、完全ＦＯＶ ＥＰＩ像を再構成する最終段階において利用され、ｍＳＥＮＳＥ解（ｍＳＥＮＳＥは、ジーメンススキャナによって実施される正則化されていないＳＥＮＳＥアルゴリズムである）と比較されている。

図１４は、４Ｄ−ＵＷＲ−ＳＥＮＳＥ法を用いてどのようにｍＳＥＮＳＥ再構成アーチファクトが取り除かれるか、軸位断、冠状断および矢状断スライス上で例示するため、２つのｐＭＲＩ再構成アルゴリズムを比較する。ｍＳＥＮＳＥ再構成画像は、感覚および認知領域（側頭葉、前頭葉および運動皮質）内の脳の中心および境界の両方に位置する大きなアーチファクトを実際に示す；図において、アーチファクトは楕円によってその概略が示される。これは結果としてＳＮ比損失をもたらし、その結果、これらの脳の領域内における賦活検出に対して劇的な効果を生じることがある。

その後、再構成パイプラインに関わらず、ＳＰＭ５ソフトウェアを用いて完全ＦＯＶ ｆＭＲＩ画像を前処理した(http://www.fil.ion.ucl.ac. uk/spm/software/spm5/）：位置ずれ補正（realignment）、運動およびスライス収集時間差の補正、空間正則化、ならびに４ｍｍ半値全幅の等方性ガウス分布核での平滑化に関する前処理である。３２チャネルヘッドコイルで収集された解剖学的Ｔ１スキャンで機能的画像の重ね合わせ（coregistration）により、ＭＮＩ空間への解剖学的正則化が実施された。ＳＰＭ５によって提供されたＴ１ ＭＮＩテンプレートにこのスキャンを正則化することによってＭＮＩ空間への正則化のためのパラメータを推定し、その後、全ての機能的画像に適用された。

被験者レベルの分析を実施するため、刺激に関連したＢＯＬＤ反応を得るべく、一般線形モデル（ＧＬＭ）を構築した。デザイン行列は、１０の実験条件に依存し、それゆえカノニカル血液動態反応関数（ＨＲＦ）でコンボリューションされたスティック関数および時間についてのその一次導関数に対応する２１の回帰子（regressor）からなり、最後の回帰子はベースラインをモデル化する。その後、同一の収集画像にこのようなＧＬＭをフィットして、ジーメンスリコンストラクタ、ＵＷＲ−ＳＥＮＳＥおよび４Ｄ−ＵＷＲ−ＳＥＮＳＥを用いて再構成した。

運動反応およびより高い認知機能（計算、言語）について対比された推定画像は、被験者およびグループレベルでのさらなる分析を受ける。これらの２つの対比は補足的なものである。なぜなら予想された賦活化が別々の脳領域にあり、それゆえ再構成アーチファクトによって別々に破損され得るからである。

より正確には、以下を研究した：
・作業記憶およびより具体的には頭頂間溝に暗算作業が作用することから、前頭および頭頂葉に誘発活動を引き起こすと考えられている聴覚計算対聴覚文（ａＣ−ａＳ）の対比。
・右運動皮質（中心前回、中前頭回）の誘発活動が期待される左クリック対右クリック（Ｌｃ−Ｒｃ）の対比。

実際に、Ｌｃ−Ｒｃ対比は、視覚または聴覚モダリティのいずれにおいても示される、２つの運動刺激を含む複合比較を規定する。それゆえ、このような比較は、運動皮質中の側性化効果を検知することを目的とする。

感覚野（視覚的、聴覚的、運動の）またはより高度な認知機能（リーディング、計算）に関連した領域を見るとき、このような理論的枠組みに対して生じた十分に異なる状況（大きい賦活化クラスター対小さい賦活化クラスター、分散した賦活化パターン対局所的な賦活化パターン、双方間活動対一方的活動）を集約するため、これらの２つの対比を選択した。

図１５は、ａＣ−ａＳ対比について、解剖学的ＭＲＩに重ね合わせられた被験者レベルのスチューデントｔマップを示す。ｍＳＥＮＳＥ、ＵＷＲ−ＳＥＮＳＥおよび４Ｄ−ＵＷＲ−ＳＥＮＳＥを用いて、それぞれＲ＝２（図の上段）およびＲ＝４（図の下段）でデータを再構成した。神経学の表現法（左は左）を採用した。十字は最大賦活化ピークを示す。

最も有意なスライスおよびＲ＝２について、全てのｐＭＲＩ再構成アルゴリズムで左頭頂葉および前頭葉の誘発活動を見出すことに成功した。しかしながら、Ｒ＝４については、ＵＷＲ−ＳＥＮＳＥおよび４Ｄ−ＵＷＲ−ＳＥＮＳＥ−そして選択的に後者−のみが、ｍＳＥＮＳＥアルゴリズムによって失われた、暗算により引き起こされた信頼性のある前頭活動を取り出すことを可能にする。定量的な観点から、提案された４Ｄ−ＵＷＲ−ＳＥＮＳＥアルゴリズムがより大きなクラスターを見出す。その極大値は、表３において説明したような、ｍＳＥＮＳＥおよびＵＷＲ−ＳＥＮＳＥを用いて得られたものより有意である。Ｒ＝２についてもっとも有意なクラスターに関して、どのような再構成アルゴリズムに対しても、そのピーク位置は安定した状態を保つ。しかしながら、それらの有意水準を調べると、ＵＷＲ−ＳＥＮＳＥをｍＳＥＮＳＥとを比較したとき、ウェーブレット基底の正則化の効果を最初に測定することができ、その後、４Ｄ−ＵＷＲ−ＳＥＮＳＥの結果を見たとき、時間的正則化および３Ｄウェーブレット分解のさらなるプラス効果を測定できる。これらの効果が、Ｒ＝４について立証もされている。表３は、（α＝０．０５の有意水準での多重比較のため修正された）ａＣ−ａＳ対比について被験者レベルでの有意な統計結果を示し、この結果は、ｐ値がαと等しいかより小さい場合に、帰無仮説が棄却されることを意味する。

統計的有意性検定の“ｐ値”を時間正則化において用いられる“形状パラメータ”、21式および22式を参照、またはＭＣＭＣアルゴリズムを例示するために用いられるパラメータ“ｐ”と混同すべきでない。

図１６は、Ｒ＝２でのａＣ−ａＳ対比のために検知された被験者間の賦活化のばらつきを例示する。実際、別のパイプライン（Ｒ＝２）で再構成された被験者レベルのスチューデントｔマップと比較すると、第２の被験者に対して期待される領域における賦活化クラスターをｍＳＥＮＳＥアルゴリズムが検出できないことが観測できる。それとは対照的に、４Ｄ−ＵＷＲ−ＳＥＮＳＥ法は、第１の被験者に対するものと正確には同一の位置にはないが、より一貫した活動を取り出す。

図１７は、Ｌｃ−Ｒｃ対比について、解剖学的ＭＲＩに重ね合わされた被験者レベルのスチューデントｔマップを示す。それぞれｍＳＥＮＳＥ、ＵＷＲ−ＳＥＮＳＥおよび４Ｄ−ＵＷＲ−ＳＥＮＳＥを用いてデータを再構成した。

全ての再構成法が右中心前回で期待される活動を取り出すことを可能にすることがわかる。しかしながら、統計的な結果をより注意深く見ると（表４を参照）、ＵＷＲ−ＳＥＮＳＥおよびより好ましくは４Ｄ−ＵＷＲ−ＳＥＮＳＥアルゴリズムは、右中前頭回のさらなるクラスターを取り出す。Ｒ＝４で収集されたデータに関し、同一のＬｃ−Ｒｃ対比が類似の活動を、すなわち同一領域において取り出す。図の最下部からわかるように、ｐＭＲＩ再構成が発明の方法により実施されたときにこの活動が高まる。

表４の定量的な結果は、図において観測できるものを数値的に裏付ける：４Ｄ−ＵＷＲ−ＳＥＮＳＥアルゴリズムを用いてＲ＝２およびＲ＝４の両方で、より高い局所的なｔスコアを有する、より大きなクラスターが検出される。より正確には、表４は（ｐ値がαより小さい、またはαに等しい場合、帰無仮説が棄却されることを意味する、α＝０．０５の有意水準での多重比較のために補正された）Ｌｃ−Ｒｃ対比についての被験者レベルでの有意な統計的結果を示す。

図１８は、このような運動の対比についての被験者間のばらつきに対する提案されたｐＭＲＩパイプラインの信頼性について説明する。感覚機能は、より大きなＢＯＬＤ効果（より高いＳＮ比）を生じ、より安定に現れることが期待されるため、Ｒ＝４で比較を行う。別のｐＭＲＩアルゴリズムを用いて再構成した２つの被験者レベルのスチューデントｔマップを比較した。第２の被験者について、右運動皮質のいかなる賦活クラスターもｍＳＥＮＳＥアルゴリズムが検知できないことを観測できる。それとは対照的に、４Ｄ−ＵＷＲ−ＳＥＮＳＥ法は、このような第２の被験者について期待される領域により一貫した活動を取り出す。

要約すると、これら２つの対比に関して、統計的有意性（クラスターの個数、クラスターの範囲、ピーク値など）の観点から、４Ｄ−ＵＷＲ−ＳＥＮＳＥアルゴリズムは常に代替再構成法よりも優れているが、信頼性の観点からも優れている。

被験者間の解剖学的および機能的なばらつきのため、集団レベルでの信頼性と再現性のある結論を導き出すためにグループレベルの分析が必要である。この検証のため、我々が以前、被験者レベルで調査した対比マップに関して１５の健康な被験者に関する変量効果分析（ＲＦＸ）が実施された。より正確には、被験者レベルの対比画像（例えば、Ｌｃ−Ｒｃ、ａＣ−ａＳ画像）について１標本のスチューデントｔ検定がＳＰＭ５を用いて実施された。

図１９はａＣ−ａＳ対比のグループレベルのスチューデントｔマップを示す。ここで、ｍＳＥＮＳＥ、ＵＷＲＳＥＮＳＥおよび４Ｄ−ＵＷＲ−ＳＥＮＳＥを用いてＲ＝２およびＲ＝４についてデータを再構成した。神経学の表現法を用いた。矢符は全体の最大賦活化ピークを示す。
これらのマップは、再構成手法に関係なく、より良いＳＮ比を考慮すると、Ｒ＝２で収集されたデータセットでより大きくてより有意な賦活化が見られることを例示する。第２に、Ｒ＝２について、目視による検査は、別々のリコンストラクタにわたるより大きなクラスター範囲および安定したクラスターについての有意性レベルの利得に加えて、４Ｄ−ＵＷＲ−ＳＥＮＳＥアルゴリズムのみが、頭頂葉皮質（軸位断のＭＩＰスライスを参照）における有意な双方間の賦活化を取り出すことを可能にすることを裏付ける。Ｒ＝４について、図の下部を見ると、類似の結論を引き出すことができる。表５において、Ｒ＝２およびＲ＝４について補足的な結果を入手でき、このような目視による比較を数値的に裏付ける：

・クラスターの範囲は強度のオーダーに反比例して減少するため、使用する再構成法がどのようなものであっても、Ｒ＝２を用いて、特にクラスターレベルで統計的性能は、はるかに有意である。
・ボクセルおよびクラスターレベルの結果は、ｍＳＥＮＳＥまたはＵＷＲ−ＳＥＮＳＥの代わりに４Ｄ−ＵＷＲ−ＳＥＮＳＥ法を用いて高められる。

表５ （ｐ＝０．０５で多重比較のために補正された）ａＣ−ａＳ対比についてグループレベルでの有意な統計的結果である。

図２０は、Ｌｃ−Ｒｃ対比に関して、Ｒ＝２およびＲ＝４についてグループレベルのＭＩＰの結果を説明する。使用する短縮ファクタＲがどのようなものであっても、運動皮質にはるかに空間的に広がった賦活化領域を我々のパイプラインが検知することを可能にすることが示される。ｍＳＥＮＳＥによって見出されたものと４Ｄ−ＵＷＲ−ＳＥＮＳＥ法を用いて検出したクラスターとを比較したとき、再びＲに関わらず、このような目視による検査が、表６において定量的に裏付けられる。最終的に４Ｄ−ＵＷＲ−ＳＥＮＳＥアルゴリズムがＵＷＲ−ＳＥＮＳＥアルゴリズムよりも優れており、このことは、提案された時空正則化スキームの利点を裏付けるものである。

表６ （ｐ＝０．０５で多重比較のために補正された）Ｌｃ−Ｒｃ対比についてグループレベルでの有意な統計的結果である。

参考文献

[Bertsekas, 1995]: Bertsekas, D. P. (1995). Nonlinear programming, Second Edition. Athena Scientific, Belmont, USA. 特に159-165頁。

[Blockら2007]: Block, K. T., Uecker, M., and Frahm, J. (2007). Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint. Magnetic Resonance in Medicine, 56 (7):1086-1098.

[Chaariら2008]: L.Chaari, J.-C. Pesquet, A. Benazza-Benyahia, P. Ciuciu, Autocalibrated Parallel MRI Reconstruction in the Wavelet Domain, in: IEEE International Symposium on Biomedical Imaging, Paris, France, 2008, pp. 756-759.

[Chaariら2009]: L. Chaari, J.-c.Pesquet, Ph. Ciuciu, and A. Benazza-Benyahia An Iterative Method for Parallel MRI SENSE-based Reconstruction in the Wavelet Domain,. arXiv:0909.0368v1 [math.OC]

[Chauxら2007]: Chaux, C., Combettes, P., Pesquet, J.-C., and Wajs, V. (2007). A variational formulation for frame-based inverse problems. Inverse Problems, 23(4):1495-1518.

[CombettesおよびPesquet, 2008]: Combettes, P. L. and Pesquet, J.-C. (2008). A proximal decomposition method for solving convex variational inverse problems. Inverse Problems, 24 (6):27.

[Daubechiesら2004]: I. Daubechies, M. Defrise, C. DeMol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Communications on Pure and Applied Mathematics 57(11) (2004) 1413-1457.

[Dempster, 1997]: Dempster, A., Laird, A., and Rubin, D. (1977). Maximum likelihood from incomplete data via the EM algorithm (with discussion). Journal of the Royal Statistical Society, Series B, 39:1-38.

[Dobigeonら2007]: Dobigeon, N. and Tourneret, J.-Y. (2007). Truncated multivariate Gaussian distribution on a simplex. Technical report, University of Toulouse.

[Donoho, 1995]: D. Donoho. De-noising by soft-thresholding. IEEE Transactions on Information Theory, 41(3):613-627, 1995.

[Geman, 1984]: S. Geman and D. Geman, “Stochastic relaxation, Gibbs distribution and the Bayesian restoration of image, " IEEE Trans. Pattern Anal. Mach. Intell., vol. 6, pp. 721-741, 1984.

[Griswoldら2002]: M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nitlka, V. Jellus, J. Wang, B. Kiefer, A. Haase, Generalized autocalibrating partially parallel acquisitions GRAPPA, Magnetic Resonance in Medicine 47 (6) (2002) 1202-1210.

[Guptaら1997]: Gupta, A. K. and Song, D. (1997). Journal of Statistical Planning and Inference. Volume 60, Pages 241-260.

[Hogeら2005]: W. S. Hoge, D. H. Brooks, B. Madore, W. E. Kyriakos, A tour of accelerated parallel MR imaging from a Iinear systems perspective, Concepts in Magnetic Resonance 27A (1) (2005) 17-37.

[Pesquet-Popescu, Pesquet] Ondeletles et applications, Techniques de l'lngenieur, traite Telecoms, TE 5 215.

[Pinelら2007]: P. Pinel, B. Thirion, S. Meriaux, A. Jobert, J. Serres， D. Le Bihan, J.-B. Poline, and S. Dehaene, “Fast reproducible identification and large-scale databasing of individual functional cognitive networks" BMC Neurosci., vol.8, no. 1, pp. 91, Oct. 2007.

[Pruessmannら1999]: K. P. Pruessmann，M .Weiger，M. B. Scheidegger, P. Boesiger, SENSE: sensitivity encoding for fast MRI, Magnetic Resonance in Medicine 42 (5) (1999) 952-962.

[Rajら2007]: Raj, A., Singh, G., Zabih, R., Kressler, B., Wang, Y., Schuff, N., and Weiner, M. (2007). Bayesian parallel imaging with edge-preserving priors. Magnetic Resonance in Medicine, 57(1):8-21.

[Robert, 1995]: “Simulation of truncated normal variables," Statistics and Computing, Vol. 5, pp. 121-125, 1995.

[Serra, 1982] Serra, J. (1982). Image Analysis and Mathematical Morphology. Academic Press, London.

[Sodicksonら1997]: Sodickson, D. K. and Manning，W.J. (1997). Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays. Magnetic Resonance in Medicine, 38(4):591-603.

[Yingら2004]: L. Ying, D. Xu, Z.-P. Liang, On Tikhonov Regularization for image reconstruction in parallel MRI, in: IEEE Engineering in Medicine and Biology Society, San Francisco, USA, 2004, pp. 1056-1059.

Claims (24)

1. − 既知のまたは推定された感度マップおよび雑音共分散行列を有するそれぞれの受信アンテナから人体の一連の基本磁気共鳴画像であって、ｋ空間でアンダーサンプリングされた基本磁気共鳴画像を収集する工程と、
   − 前記人体の磁気共鳴画像の正則化再構成を行う工程とを備え、
   前記磁気共鳴画像の正則化再構成を行う工程は、
   − 前記基本磁気共鳴画像に基づき再構成された再構成画像の尤度を表す誤差項、
   − 再構成画像の全変分ノルムである少なくとも１つの空間領域ペナルティ項および
   − 前記再構成画像のフレーム係数の実際の統計的分布と前記係数の事前分布との間の偏差を表すフレームペナルティ項
   を備えた評価関数を最小化することによって離散フレーム空間で実施され、
   再構成画像のフレーム係数の前記事前分布は、前記人体の補助的な磁気共鳴画像に基づいて推定される、人体のパラレル磁気共鳴撮像法。
2. 前記誤差項は、前記基本磁気共鳴画像に基づき再構成された前記再構成画像の負の対数尤度を表す、請求項１に記載の磁気共鳴撮像法。
3. 前記人体の磁気共鳴画像の正則化再構成を行う前記工程は、前記基本磁気共鳴画像とフレーム係数の前記事前分布を考慮して、人体の画像を規定する一連のフレーム係数の完全な事後分布を、前記フレーム空間において、最大化することによって行われる、請求項１または２に記載の磁気共鳴撮像法。
4. 前記人体の前記補助的な磁気共鳴画像は、前記基本磁気共鳴画像から再構成される、請求項１ないし３のいずれか１つに記載の磁気共鳴撮像法。
5. 前記人体の前記補助的な磁気共鳴画像は、感度エンコーディング−SENSE−アルゴリズ
   ムを用いて前記基本磁気共鳴画像から再構成される、請求項４に記載の磁気共鳴撮像法。
6. 前記人体の前記補助的な磁気共鳴画像は、
   − 非正則化SENSEアルゴリズム、
   − 画像空間で正則化されたSENSEアルゴリズムおよび
   − k空間において正則化されたSENSEアルゴリズム
   から選択されたアルゴリズムを用いた前記基本磁気共鳴画像から再構成される、請求項５に記載の磁気共鳴撮像法。
7. 前記フレーム係数の汎用ガウス−ラプラス事前統計的分布を想定し、そして前記分布のパラメータは、最尤推定量または事後平均推定量を用いて、前記人体の前記補助的な磁気共鳴画像に基づき推定される、請求項１ないし６のいずれか１つに記載の磁気共鳴撮像法。
8. 前記誤差項は、二次平均誤差項である、請求項１ないし７のいずれか１つに記載の磁気共鳴撮像法。
9. 前記基本磁気共鳴画像は三次元画像であり、そして磁気共鳴画像の正則化再構成を行う前記工程は、離散三次元フレーム空間で実施される、請求項１ないし８のいずれか１つに記載の磁気共鳴撮像法。
10. 前記基本磁気共鳴画像は、撮像すべき被験者をスライスした二次元基本磁気共鳴画像を積層させることによって得られる請求項９に記載の磁気共鳴撮像法。
11. 磁気共鳴画像の正則化再構成を行う前記工程は、冗長ウェーブレットフレーム表現に基づく、請求項１ないし１０のいずれか１つに記載の磁気共鳴撮像法。
12. 磁気共鳴の正則化再構成を行う前記工程は、非冗長ウェーブレット表現に基づく、請求項１ないし１０のいずれか１つに記載の磁気共鳴撮像法。
13. − 既知または推定された感度マップおよび雑音共分散行列を有するそれぞれの受信アンテナから人体の一連の時系列の基本磁気共鳴画像であって、ｋ空間でアンダーサンプリングされた基本磁気共鳴画像を収集する工程と、
    − 前記人体の一時系列の磁気共鳴画像の正則化再構成を行う工程とを備え、
    一時系列の基本磁気共鳴画像の正則化再構成を行う工程は、
    − 対応する前記基本磁気共鳴画像に基づき再構成された各再構成画像の尤度を表す誤差項および
    − 再構成画像の全変分ノルムである少なくとも１つの空間領域ペナルティ項、
    − 前記再構成画像のフレーム係数の実際の統計的分布と前記係数の事前分布との間の偏差を表す、フレームペナルティ項および
    − その系列の連続した画像間の画素ごとまたはボクセルごとの差を表す、エッジ保存時間的ペナルティ項
    を備えた微分不可能な評価関数を最小化することによって実施され、
    再構成画像のフレーム係数の前記事前分布は、前記人体の補助的な磁気共鳴画像に基づき推定される、人体のダイナミックおよびパラレル磁気共鳴撮像法の実施方法。
14. 前記時間的ペナルティ項が、凸エッジ保存関数に基づく、請求項１３に記載の方法。
15. 前記時間的ペナルティ項が、p≧１のＬpノルムに基づく、請求項１４に記載の方法。
16. 前記時間的ペナルティ項は、第１部分時間的ペナルティ項および第２部分時間的ペナルティ項の和によって与えられ、
    − 第１部分時間的ペナルティ項は、その系列の各偶数画像と先行する奇数画像との間の
    ピクセルごとまたはボクセルごとの差を表し、
    − 第２部分時間的ペナルティ項は、その系列の各奇数画像と先行する偶数画像との間のピクセルごとまたはボクセルごとの差を表し、
    前記評価関数は、前記第１および第２部分時間的ペナルティ項に対し、近接演算子を用いることによって最小化される、請求項１３ないし１５のいずれか１つに記載の方法。
17. 前記基本磁気共鳴画像は三次元画像であり、前記離散フレーム空間は、離散三次元フレーム空間である、請求項１６に記載の方法。
18. 最尤推定量を用いて前記時間的ペナルティ項の重み付けパラメータを自動的に決定する工程をさらに備えた、請求項１３ないし１７のいずれか１つに記載の方法。
19. 再構成すべき画像の、各ピクセル若しくはボクセルまたは一連の隣接したピクセル若しくはボクセルに対する時間的ペナルティ項の前記重み付けパラメータを推定する工程を備えた、請求項１８に記載の方法。
20. 前記時間的ペナルティ項がＬpノルムに基づき、時間的ペナルティ項の前記重み付けパ
    ラメータおよびパラメータｐが前記最尤推定量を用いて共に決定される、請求項１８または１９に記載の方法。
21. 時間的ペナルティ項の前記重み付けパラメータおよびパラメータｐは、制約ｐ≧１の下で前記最尤推定量を用いて共に決定される、請求項２０に記載の方法。
22. 前記誤差項は、基準とみなされる基本磁気共鳴画像に対する前記基本磁気共鳴画像のそれぞれの剛体変換を規定する幾何学的パラメータに依存し、そして一時系列の基本磁気共鳴画像の正則化再構成を行う前記工程は、前記幾何学的パラメータについて前記関数を最小化することによっても実施される、請求項８ないし２０のいずれか１つに記載の方法。
23. 前記基本磁気共鳴画像は、エコープラナー撮像法によって収集される、請求項１ないし２２のいずれか１つに記載の方法。
24. 前記基本磁気共鳴画像は、４より大きいかまたは４に等しい減少係数でアンダーサンプリングされる、請求項１ないし２３のいずれか１つに記載の方法。

Images