本発明の実施形態に係る磁気共鳴イメージング装置について添付図面を参照して説明する。
図1は本発明の実施形態に係る磁気共鳴イメージング装置の構成図である。
磁気共鳴イメージング装置20は、静磁場を形成する筒状の静磁場用磁石21、この静磁場用磁石21の内部に設けられたシムコイル22、傾斜磁場コイル23及びRFコイル24を備えている。
また、磁気共鳴イメージング装置20には、制御系25が備えられる。制御系25は、静磁場電源26、傾斜磁場電源27、シムコイル電源28、送信器29、受信器30、シーケンスコントローラ31及びコンピュータ32を具備している。制御系25の傾斜磁場電源27は、X軸傾斜磁場電源27x、Y軸傾斜磁場電源27y及びZ軸傾斜磁場電源27zで構成される。また、コンピュータ32には、入力装置33、表示装置34、演算装置35及び記憶装置36が備えられる。
静磁場用磁石21は静磁場電源26と接続され、静磁場電源26から供給された電流により撮像領域に静磁場を形成させる機能を有する。尚、静磁場用磁石21は超伝導コイルで構成される場合が多く、励磁の際に静磁場電源26と接続されて電流が供給されるが、一旦励磁された後は非接続状態とされるのが一般的である。また、静磁場用磁石21を永久磁石で構成し、静磁場電源26が設けられない場合もある。
また、静磁場用磁石21の内側には、同軸上に筒状のシムコイル22が設けられる。シムコイル22はシムコイル電源28と接続され、シムコイル電源28からシムコイル22に電流が供給されて静磁場が均一化されるように構成される。
傾斜磁場コイル23は、X軸傾斜磁場コイル23x、Y軸傾斜磁場コイル23y及びZ軸傾斜磁場コイル23zで構成され、静磁場用磁石21の内部において筒状に形成される。傾斜磁場コイル23の内側には寝台37が設けられて撮像領域とされ、寝台37には被検体Pがセットされる。RFコイル24にはガントリに内蔵されたRF信号の送受信用の全身用コイル(WBC: whole body coil)や寝台37や被検体P近傍に設けられるRF信号の受信用の局所コイルなどがある。
また、傾斜磁場コイル23は、傾斜磁場電源27と接続される。傾斜磁場コイル23のX軸傾斜磁場コイル23x、Y軸傾斜磁場コイル23y及びZ軸傾斜磁場コイル23zはそれぞれ、傾斜磁場電源27のX軸傾斜磁場電源27x、Y軸傾斜磁場電源27y及びZ軸傾斜磁場電源27zと接続される。
そして、X軸傾斜磁場電源27x、Y軸傾斜磁場電源27y及びZ軸傾斜磁場電源27zからそれぞれX軸傾斜磁場コイル23x、Y軸傾斜磁場コイル23y及びZ軸傾斜磁場コイル23zに供給された電流により、撮像領域にそれぞれX軸方向の傾斜磁場Gx、Y軸方向の傾斜磁場Gy、Z軸方向の傾斜磁場Gzを形成することができるように構成される。
RFコイル24は、送信器29及び受信器30の少なくとも一方と接続される。送信用のRFコイル24は、送信器29からRF信号を受けて被検体Pに送信する機能を有し、受信用のRFコイル24は、被検体P内部の原子核スピンのRF信号による励起に伴って発生したNMR信号を受信して受信器30に与える機能を有する。
図2は図1に示すRFコイル24の詳細構成の一例を示す図である。
図2に示すようにRFコイル24は、筒状の全身用(WB: whole-body)コイル24aとフェーズドアレイコイル24bを備えている。フェーズドアレイコイル24bは、複数のコイル要素24cを備えており、被検体Pの体表側と背面側とにそれぞれ複数のコイル要素24cが配置される。
一方、受信器30は、デュプレクサ30a,アンプ30b、切換合成器30cおよび受信系回路30dを備えている。デュプレクサ30aは、送信器29、WBコイル24aおよびWBコイル24a用のアンプ30bと接続される。アンプ30bは、各コイル要素24cおよびWBコイル24aの数だけ設けられ、それぞれ個別に各コイル要素24cおよびWBコイル24aと接続される。切換合成器30cは、単一または複数個設けられ、切換合成器30cの入力側は、複数のアンプ30bを介して複数のコイル要素24またはWBコイル24aと接続される。受信系回路30dは、各コイル要素24cおよびWBコイル24aの数以下となるように所望の数だけ設けられ、切換合成器30cの出力側に設けられる。
WBコイル24aは、RF信号の送信用のコイルとして用いることができる。また、NMR信号の受信用のコイルとして各コイル要素24cを用いることができる。さらに、WBコイル24aを受信用のコイルとして用いることもできる。
このため、デュプレクサ30aは、送信器29から出力された送信用のRF信号をWBコイル24aに与える一方、WBコイル24aにおいて受信されたNMR信号を受信器30内のアンプ24dを経由して切換合成器30cに与えるように構成されている。また、各コイル要素24cにおいて受信されたNMR信号もそれぞれ対応するアンプ24dを経由して切換合成器30cに出力されるように構成されている。
切換合成器30cは、コイル要素24cやWBコイル24aから受けたNMR信号の合成処理および切換を行って、対応する受信系回路30dに出力するように構成されている。換言すれば、受信系回路30dの数に合わせてコイル要素24cやWBコイル24aから受けたNMR信号の合成処理および切換が切換合成器30cにおいて行われ、所望の複数のコイル要素24cを用いて撮影部位に応じた感度分布を形成して様々な撮影部位からのNMR信号を受信できるように構成されている。
ただし、コイル要素24cを設けずに、WBコイル24aのみでNMR信号を受信するようにしてもよい。また、切換合成器30cを設けずに、コイル要素24cやWBコイル24aにおいて受信されたNMR信号を直接受信系回路30dに出力するようにしてもよい。さらに、より多くのコイル要素24cを広範囲に亘って配置することもできる。
一方、制御系25のシーケンスコントローラ31は、傾斜磁場電源27、送信器29及び受信器30と接続される。シーケンスコントローラ31は傾斜磁場電源27、送信器29及び受信器30を駆動させるために必要な制御情報、例えば傾斜磁場電源27に印加すべきパルス電流の強度や印加時間、印加タイミング等の動作制御情報を記述したシーケンス情報を記憶する機能と、記憶した所定のシーケンスに従って傾斜磁場電源27、送信器29及び受信器30を駆動させることによりX軸傾斜磁場Gx、Y軸傾斜磁場Gy,Z軸傾斜磁場Gz及びRF信号を発生させる機能を有する。
また、シーケンスコントローラ31は、受信器30におけるNMR信号の検波及びA/D (analog to digital)変換により得られた複素データである生データ(raw data)を受けてコンピュータ32に与えるように構成される。
このため、送信器29には、シーケンスコントローラ31から受けた制御情報に基づいてRF信号をRFコイル24に与える機能が備えられる一方、受信器30には、RFコイル24から受けたNMR信号を検波して所要の信号処理を実行するとともにA/D変換することにより、デジタル化された複素データである生データを生成する機能と生成した生データをシーケンスコントローラ31に与える機能とが備えられる。
さらに、磁気共鳴イメージング装置20には、被検体PのECG (electro cardiogram)信号を取得するECGユニット38が備えられる。ECGユニット38により取得されたECG信号はシーケンスコントローラ31を介してコンピュータ32に出力されるように構成される。尚、拍動を心拍情報として表すECG信号の代わりに拍動を脈波情報として表す脈波同期(PPG: peripheral pulse gating)信号を取得することもできる。
また、コンピュータ32の記憶装置36に保存されたプログラムを演算装置35で実行することにより、コンピュータ32には各種機能が備えられる。ただし、プログラムの少なくとも一部に代えて、各種機能を有する特定の回路を磁気共鳴イメージング装置20に設けてもよい。
図3は、図1に示すコンピュータ32の機能ブロック図である。
コンピュータ32の演算装置35は、記憶装置36に保存されたプログラムを実行することにより撮像条件設定部40、データ処理部41及びAFI処理条件設定部42として機能する。また、記憶装置36は、k空間データベース43、画像データベース44、位相データベース45及びAFI条件データベース46として機能する。データ処理部41は、PI (parallel imaging)処理部41A及びAFI処理部41Bを有する。
撮像条件設定部40は、入力装置33からの指示情報に基づいてパルスシーケンスを含む撮像条件を設定し、設定した撮像条件をシーケンスコントローラ31に出力することによってシーケンスコントローラ31を駆動制御させる機能を有する。
特に、撮像条件設定部40は、k-spaceにおいて波数方向に非対称となるようにMRデータをサンプンリングするAFIの撮像条件及び複数のコイル要素24cを用いてデータを受信するPIの撮像条件を設定できるように構成されている。AFI及びPIは、磁気共鳴血管撮影法(MRA: magnetic resonance angiography)や拡散強調イメージング(DWI: diffusion weighted imaging)を含む様々なイメージングに適用することができる。
尚、PIには、SENSE (Sensitivity Encoding)、SMASH (Shimadzu minimum angle shot)、GRAPPA (generalized autocalibrating partially parallel acquisition)等の種類がある。SENSEではr-spaceにおいてデータ処理が実行されるのに対してSMASHではk-spaceにおいてデータが処理される。GRAPPAは、SMASHを発展させたもので、PI特有の後処理である展開処理及び画像再構成処理をコイル要素24c毎のk-spaceデータに実行した後に展開処理後のコイル要素24c毎の画像データを合成する手法である。
AFIを行うための撮像条件としては、2次元(2D: two dimensional)のサンプリングであれば、k-spaceデータの読み出し(readout)方向及び位相エンコード(phase encode)方向のうち一方向kの波数方向におけるサンプリング領域(-Kc≦k≦Kmax)がある。サンプリング領域の境界Kcは、例えば、AFI処理条件設定部42からの設定情報に従って可変設定することができる。非サンプリング領域は、k方向の正側及び負側のいずれでもよいが、ここでは負側を非サンプリング領域とする場合を例に説明する。
或いは、サンプリング領域の境界Kcは、プレスキャンによって収集されたデータに基づいて決定することもできる。この場合、プレスキャン用のデータ収集条件は撮像条件設定部40において設定される。具体的には、サンプリング領域の境界Kcを徐々に変化させてデータ収集を行う条件がプレスキャン用のデータ収集条件として設定される。
図4は、AFIにおけるデータサンプリング領域の境界Kcをプレスキャンによって決定する方法を説明する図である。
図4(A)に示すように、データサンプリング領域の境界KcをKc1, Kc2, Kc3,...と徐々に変化させてデータ収集を行う条件がプレスキャン用に設定される。この境界Kcを求めるためのプレスキャンをKc-PREPスキャンと呼ぶことにすると、Kc-PREPスキャンのデータ収集条件は、イメージングスキャンのデータ収集条件となるべく一致させることが、より適切な境界Kcを求めることに繋がる。一方、データ収集時間の短縮化のために位相エンコードをゼロにして、つまり1次元(1D: one-dimensional)の投影データとしてデータサンプリング領域の境界Kcを変化させてもよい。
Kc-PREPスキャンにより収集された非対称サンプリングデータに対して後述するAFI処理を行うと、図4(B)に示すように、互いに異なるデータサンプリング領域の境界Kc1, Kc2, Kc3,...にそれぞれ対応する2D画像I (Kc1), I (Kc2), I (Kc3), ...が得られる。そうすると、ユーザから入力装置33を通じた画像選択情報或いは閾値処理等の画像処理により、適切な画質の画像I (Kc_opt)を選択することができる。そして、この画像データの選択情報に基づいて、適切な画質の画像I (Kc_opt)に対応するデータサンプリング領域の境界Kc_optをイメージングスキャン用の撮像条件として決定することができる。
そして、図4(C)に示すように、選択されたデータサンプリング領域の境界Kc_optを撮像条件としてイメージングスキャンが実行される。
更に、撮像条件設定部40は、この他のイメージングスキャン用の撮像条件の設定やイメージングスキャンで収集されたデータの処理に必要なデータを取得するためのプレスキャン用のデータ収集条件を設定する機能を備えている。例えば、静磁場のシミングを行う場合には、シミング用のプレスキャンのデータ収集条件が設定される。シミング用のプレスキャンで収集されたMRデータからは、シムコイル電源28からシムコイル22に供給すべき電流を決定するための磁場分布マップが測定される。別の例としては、PIによって収集されたデータに対する展開処理に必要なコイル要素24cの感度マップを測定するためのプレスキャンがある。
データ処理部41は、シーケンスコントローラ31からAFI及びPI用の条件によって収集された生データを受けてk空間データベース42に形成されたk-spaceにk-spaceデータとして配置する機能と、k-spaceデータに対するAFI処理及びPIに起因するPI処理を実行することによって診断用の画像データを生成する機能とを有する。AFI処理は、k-spaceで波数方向に非対称にサンプンリングされたk-spaceに基づいて、位相補正処理を伴って対称にサンプリングされたデータから生成される画像データと同等な画像データを生成する処理である。PI処理及びAFI処理によって生成された診断用の画像データは、表示装置34に表示させることができる。
データ処理部41のPI処理部41Aは、PI処理を実行する機能を有し、AFI処理部41BはAFI処理を実行する機能を有する。PI処理及びAFI処理では、k-spaceデータとr-spaceデータ間の変換が必要となる場合がある。そこで、PI処理部41A及びAFI処理部41Bには、k-spaceデータとr-spaceデータ間の変換処理機能が設けられるとともに、k空間データベース42からのk-spaceデータの読み込み、k空間データベース42へのk-spaceデータの書き込み、画像データベース43からのr-spaceデータの読み込み及び画像データベース43へのr-spaceデータの書き込みを行う機能が備えられる。
このため、k空間データベース42には、k-spaceデータが保存される一方、画像データベース43にはr-spaceデータが保存される。
PI処理には、折り返しの発生した画像データを展開することによって折り返しを除去する展開処理及びコイル要素24c毎の画像データを合成する合成処理がある。一方、AFI処理には、複素共役対称性を利用したデータ補填のためのk-spaceデータのフィルタ処理、r-spaceデータの位相補正処理、位相補正エラーの低減処理等の処理が挙げられる。
図5は、図3に示すデータ処理部41において実行される第1のAFI処理の例を示す図である。
図5は、k-spaceのkx方向に非対称なk-spaceデータがサンプリングされた例を示しているがky方向に非対称なk-spaceデータがサンプリングされた場合も同様である。また、3次元(3D: three-dimensional)k-spaceデータを2D方向に非対称にサンプリングしてAFI処理を行うことも可能である。以降では、説明簡易化のため、k方向に1Dの非対称サンプリングを行う場合を例に説明する。
尚、データを非対称にサンプリングする方向を読み出し(read-out)方向とすればエコー時間(TE: echo time)の短縮が可能である。一方、データを非対称にサンプリングする方向をエンコード方向とすれば撮像時間の短縮が可能である。
図5に示すように、AFI法では、k-spaceに中心に関して非対称となるようにMRデータが収集され、k-spaceにk-spaceデータとして配置される。このため、k-spaceには、k=0に関して非対称なサンプリング領域が形成され、部分的なk-spaceデータが非対称サンプリング領域に充填される。図5は1Dで、一方の高周波側のデータが欠落した横軸方向の-Kc≦k≦Kmaxの範囲におけるk-spaceデータが、部分的原k-spaceデータSorig(k)として配置された例を示している。
第1のAFI処理では、まず式(1)に示すように部分的原k-spaceデータSorig(k) (-Kc≦k≦Kmax)はhomodyne filter Hhomo(k)が掛けられた後、FTによって原r-spaceデータIhomo(x)が生成される。
homodyne filterは、k-spaceの無データ部分に複素共役データを充填するのと等価なフィルタである。すなわち、homodyne filterによる部分的原k-spaceデータSorig(k)のwindowingによって、複素共役対称を利用したk-spaceデータの補填処理が実行される。homodyne filterは、k-spaceデータの対称サンプリング部分と非対称サンプリング部分とで重みが異なるwindow関数となる。
[数1]
Ihomo(x) = FT{Hhomo(k)Sorig(k)} (1)
一方、式(2)に示すように部分的原k-spaceデータSorig(k)にはLPF (low pass filter) Hlow(k)が掛けられ、k-spaceにおいて対称な低周波領域の原k-spaceデータSlow(k) (-Kc≦k≦Kc)が抽出される。そして、抽出された低周波領域の原k-spaceデータのFTによって、低周波データに対応するlow-pass filtered r-spaceデータIlow(x)が生成される。
[数2]
Ilow(x) = FT{Hlow(k)Sorig(k)} (2)
そして、式(3-1)及び式(3-2)に示すようにlow-pass filtered r-spaceデータIlow(x)から低周波領域における位相分布Φlow(x)を求め、低周波領域における位相分布Φlow(x)を用いて原r-spaceデータIhomo(x)の位相補正が実行された後、虚数部を除去するREAL化処理が行われる。これにより、AFI処理後のr-spaceデータIcor(x)が生成される。
[数3]
exp{-iΦlow(x)} = Ilow *(x)/abs{Ilow(x)} (3-1)
Icor(x) = Re[Ihomo(x)exp{-iΦlow(x)}] (3-2)
但し、式(3-2)においてRe[]は、複素数の実部を出力する関数である。
すなわち、AFIは、複素共役の対称性を利用して非サンプリング部分のk-spaceデータを埋めるイメージング法であるが、実際には、位相分布が無視できず、虚部がゼロとならないため対称性が成立しない。そこで、位相分布に基づく位相補正によって複素共役の対称性の確保が図られる。しかしながら、この位相補正によっても補正エラーが発生する場合がある。
そこで、位相補正後のr-spaceデータIcor(x)に対し位相補正のエラーを低減させるためのループ処理が所定回数実行される。このループ処理は、位相補正にエラー(誤差)がなければ、k-spaceでの複素共役が対称となるため位相補正後のr-spaceデータIcor(x)の虚部がゼロになるという原理に基づいている。従って、ループ処理は、r-spaceデータIcor(x)の虚部をゼロに収束させる収束計算となる。そして、ループ処理を繰り返す収束計算によってhomodyne filter処理に起因して生じる位相補正のエラーが低減される。
具体的には、式(4)に示すように位相補正後の実部のみのr-spaceデータIcor(x)をループ処理における初期のr-spaceデータIj(x)とし、r-spaceデータIj(x)の位相を位相補正前の位相に戻す処理が実行される。これにより位相を戻したr-spaceデータI'j(x)が生成される。
[数4]
I'j(x) = Ij(x)exp{iΦlow(x)} (4)
次に、式(5)に示すように位相を戻したr-spaceデータI'j(x)がIFTによりk-spaceデータSj(k)に変換される。
[数5]
Sj(k) = IFT{I'j(x)} (5)
次に、式(6)に示すようにk-spaceデータSj(k)のうち-Kc≦k≦Kmaxの範囲が部分的原k-spaceデータSorig(k)からデータ切出しフィルタHmerge(k)によって切出されたk-spaceデータと置換される。換言すれば、式(5)で得られたk-spaceデータSj(k)のうちのk<-Kcの範囲のk-spaceデータと部分的原k-spaceデータSorig(k)のうちの-Kc≦k≦Kmaxの範囲のk-spaceデータが合成される。
[数6]
Sj+1(k) = Hmerge(k)Sorig(k)+{1-Hmerge(k)}Sj(k) (6)
次に、式(5)で得られた合成後のk-spaceデータSj+1(k)のFTによって得られるr-spaceデータに対して位相補正が行われた後、虚数部を除去するREAL化処理が行われる。これにより、式(7)に示すように再度位相補正後の実部のみのr-spaceデータIj+1(x)が生成される。
[数7]
Ij+1(x) = Re[FT{Sj+1(k)}exp{-iΦlow(x)}] (7)
式(4)から式(7)までのループ処理は、所定の収束条件が満たされるまでj=0, 1, 2, 3, ...として繰り返される。収束条件は、例えば式(8)のように、j番目のr-spaceデータIj(x)とj+1番目のr-spaceデータIj+1(x)との差が閾値Eth未満となる場合とすることができる。
[数8]
|Ij+1(x)-Ij(x)|<Eth (8)
そして、収束条件が満たされるまでループ処理が繰り返されることによって、診断画像データとしてのr-spaceデータI(x)が生成される。このため、第1のAFI処理では、ループ処理の繰り返し回数をj回とすると、FTの回数は2+j回、IFTの回数はj回となる。このため、第1のAFI処理におけるデータ処理時間は、2+j回のFT及びj回のIFTに要する時間を含む処理時間となる。
図6は、図3に示すデータ処理部41において実行される第2のAFI処理の例を示す図である。
図6に示す第2のAFI処理では、主に部分的原k-spaceデータSorig(k)にhomodyne filter Hhomo(k)を掛ける前に位相補正処理を行う点が第1のAFI処理と異なる。このため、第1のAFI処理と同様な処理については説明を省略する。
第2のAFI処理でも、第1のAFI処理と同様に式(2)により部分的原k-spaceデータSorig(k)から低周波データに対応するlow-pass filtered r-spaceデータIlow(x)が生成される。
一方、式(9)に示すように部分的原k-spaceデータSorig(k)から全周波数領域のデータを切出すフィルタHwholeにより、非対称な原k-spaceデータが切出され、FTにより原r-spaceデータIwhole(x)が生成される。
[数9]
Iwhole(x) = FT{HwholeSorig(k)} (9)
次に、式(10)に示すようにlow-pass filtered r-spaceデータIlow(x)から求められた低周波領域における位相分布Φlow(x)を用いて原r-spaceデータIwhole(x)の位相補正が実行される。これにより位相補正後のr-spaceデータIpcor(x)が得られる。
[数10]
Ipcor(x) = Iwhole(x)exp{-iΦlow(x)}] (10)
次に、式(11)に示すように位相補正後のr-spaceデータIpcor(x)がIFTにより位相補正後のk-spaceデータSPcor(k)に変換される。
[数11]
SPcor(k) = IFT{Ipcor(x)} (11)
次に、式(12)に示すように位相補正後のk-spaceデータSPcor(k)にhomodyne filter Hhomo(k)が掛けられた後、FTによってr-spaceデータに変換され、更にr-spaceデータの虚数部を除去するREAL化処理が行われる。これにより、homodyne filter処理及び位相補正後のr-spaceデータIcor(x)が生成される。
[数12]
Icor(x) = Re[FT{Hhomo(k)SPcor(k)}] (12)
このように生成された位相補正後のr-spaceデータIcor(x)に対して第1のAFI処理と同様にループ処理が所定回数実行される。このため、第2のAFI処理では、ループ処理の繰り返し回数をj回とすると、FTの回数は3+j回、IFTの回数は1+j回となる。従って、第2のAFI処理におけるデータ処理時間は、FT及びIFTそれぞれ1回分に相当する時間だけ第1のAFI処理より長くなる。ただし、部分的原k-spaceデータSorig(k)にhomodyne filter Hhomo(k)を掛ける前に位相補正処理が実行される。このため、第2のAFI処理では、第1のAFI処理に比べてhomodyne filter処理に起因する位相補正のエラーを低減させることができる。
図7は、図3に示すデータ処理部41において実行される第3のAFI処理の例を示す図である。
図7に示す第3のAFI処理では、非対称にサンプリングされた部分的原k-spaceデータSorig(k)全体に基づいて全周波数領域の位相分布Φwhole(x)を推定し、全周波数領域における位相分布Φwhole(x)を用いて位相補正が実行される点が第2のAFI処理と異なる。すなわち、対称にサンプリングされた部分のみならず、非対称にサンプリングされた部分の部分的原k-spaceデータSorig(k)も用いることによって、全周波数領域の位相分布Φwhole(x)が位相補正用に推定される。他の点については第2のAFI処理と同様であるため詳細な説明を省略する。
第3のAFI処理では、全周波数領域における位相分布Φwhole(x)が位相補正用に用いられる。このため、第3のAFI処理では、第2のAFI処理よりも更に位相補正のエラーを低減させることができる。尚、第3のAFI処理では、FTは3+j回実行され、IFTは1+j回実行される。
上述のAFI処理方法では、位相補正用の位相分布が補正対象となるk-spaceデータ自身から推定されるが、診断画像データの生成に用いられるMRデータ以外の他のMRデータから位相補正用の位相分布を事前に求めておくこともできる。特にAFI法を用いずにデータ収集を行うと、診断画像データの生成用に非対称にサンプリングされたk-spaceデータよりも広い周波数領域におけるk-spaceデータが得られる。例えば、全周波数領域のk-spaceデータを収集すれば、全周波数領域の位相分布Φwhole(x)を求めることができる。
位相分布を自己データではなく他のデータから求める場合には、同一の被検体の同一の撮像部位から収集したデータから位相分布を求めることが精度上実用的である。すなわち、より本来の位相分布に近い位相分布を推定することができる。この場合、他のデータは、少なくとも自己データよりも高周波領域までサンプリングされていることが重要である。
例えば、静磁場のシミング用の磁場マップを測定するためのシミングシーケンスにより、ほぼfull samplingされたk-spaceデータが得られる。また、PI処理用のコイル要素24cの感度分布を測定するためにも、自己データより広いk空間(周波数領域)のデータがプレスキャンによって収集される。このため、シミングシーケンスや感度マップ取得用シーケンスにより収集されたデータを利用して全周波数領域における位相補正用の位相分布を求めることができる。
そこで、AFI処理部41Bには、プレスキャン等の所望のスキャンによって収集された他のデータに基づいて位相補正用の位相分布を求める機能、求めた位相分布を位相データベース45に書き込む機能及び位相データベース45から位相分布を読み込む機能が備えられる。このため、位相データベース45には、事前に取得されたAFI処理における位相補正用の位相分布が保存される。
位相分布を他のデータから求める場合、他のデータ及び自己データの双方を用いてもよい。例えば、シーケンスに応じた他のデータの位相シフトと、位相補正対象となるデータの位相シフトの類似度を表す指標が閾値以上であれば、他のデータから求めた全周波数領域における位相分布を位相補正用の位相分布とする。一方、位相シフト間の類似度を表す指標が閾値未満であれば、ほぼfull samplingされた他のk-spaceデータのうち、位相補正対象となる自己データにおいてサンプリングされない部分のデータをk-spaceにおけるwindowingにより抽出する。そして、抽出したk-spaceデータをr-spaceデータに変換し、得られたr-spaceデータに基づいて算出した位相分布を、自己データに基づいて算出した位相分布と合成する。
これにより、他のデータ及び自己データの双方に基づく全周波数領域における位相分布マップを求めることができる。このような位相分布の算出方法により、厳密には本来の位相分布と異なる位相分布が位相補正に用いられるものの、位相補正エラーの発生を低減させることができる。
尚、上述の例では、位相補正がk-spaceデータに対して実行されているがr-spaceデータに対して位相補正を実行してもよい。
ところで、PI処理及びAFI処理は、処理順序を最適化することによって処理時間を短縮させることができる。そのため、PI処理部41A及びAFI処理部41Bは、それぞれ処理中又は処理後の中間データを他方に与えるように構成される。例えば、AFI処理におけるhomodyne filter処理を複数のコイル要素24cに対応する各k-spaceデータに実行した後に、PI処理における展開処理及び画像データ合成処理を行うようにすれば、位相補正のエラーを低減させるためのループ処理の対象が単一のr-spaceデータとなる。換言すれば、PI処理における展開処理及び画像データ合成処理後にループ処理を行うようにすれば、ループ処理の対象がコイル要素24c間で合成された画像データのみとなる。このため、ループ処理の回数を減らし、全体のデータ処理量及びデータ処理時間を低減させることができる。
AFI処理条件設定部42は、位相補正処理に用いられる位相分布又は位相分布に影響のある撮像条件に応じてAFI処理の条件を適切な条件に設定する機能を有する。AFI処理の条件としては、上述のようにAFI処理の種類、ループ処理の回数j、サンプリング領域の境界Kcの他、homodyne filterの強度及び形状などがある。homodyne filterの強度は、非対称にサンプリングされたMR信号のゲインとなる。
AFI処理の種類は、例えば3つの変更方法をパラメータ化することによって特定することができる。第1の変更方法は、位相補正処理をhomodyne filter処理前に行うかhomodyne filter処理後に行うかという変更方法である。第2の変更方法は、データ処理範囲に対応する全周波数領域の位相分布Φwhole(x)を位相補正に用いるか非対称サンプリング領域のうち対称な部分に対応する低周波領域の位相分布Φlow(x)を位相補正に用いるかという変更方法である。第3の変更方法は、位相補正に用いる位相分布Φ(x)を補正対象となるk-spaceデータ自身、つまり自己データのみから推定するか、他のデータのみから推定するか、自己データ及び他のデータの双方から推定するかという変更方法である。
AFI処理の条件となる各種パラメータは、入力装置33の操作によってAFI処理条件設定部42に設定情報を入力することにより位相補正に用いられる位相分布又は位相分布に影響のある撮像条件に応じて可変設定することができる。また、位相分布又は撮像条件に応じた複数のAFI処理の条件を候補として予め準備しておき、入力装置33の操作によって選択できるようにすることもできる。更に、様々な撮像条件に対応する複数の撮像プロトコルにそれぞれ択一的で適切なAFI処理条件のパラメータを関連付けて保存しておくこともできる。
AFI条件データベース46には、予めシミュレーションにより、又は経験的に決定した複数の撮像条件に対応するAFI処理条件又は撮像条件に応じたAFI処理条件の候補が保存される。そして、AFI処理条件設定部42は、撮像条件設定部40から撮像プロトコルの選択情報等の撮像条件情報を取得した場合に、対応するAFI処理条件又はAFI処理条件の候補をAFI条件データベース46から取得するように構成される。この場合、例えば入力装置33の操作によって撮像プロトコルを選択すると、対応するAFI処理条件のパラメータがユーザに無意識に自動設定されるため、ユーザの利便性向上に繋がる。
但し、サンプリング領域の境界Kcについては上述したように、Kc-PREPスキャンの実行によって決定することもできる。
位相分布に影響のある撮像条件としては、イメージングシーケンスの種類、TE、撮像部位等の条件が挙げられる。例えば、FASE (fast advanced spin echo又はfast asymmetric spin echo)等のFE (field echo)系のシーケンス、TEの長いGE (gradient echo)系のシーケンス、EPI (echo planar imaging)シーケンス、定常自由歳差運動(SSFP: steady state free precession)シーケンスを用いてイメージングを行う場合には、位相補正の誤差によってアーチファクトが顕著に現れる場合がある。
また、撮像部位が血管であるMRA、DWIを行う場合や磁気共鳴イメージング装置20が高磁場装置である場合には、高精度な位相補正が要求される。MRAには、FS-BB (flow-sensitive black-blood)法やTOF (time of flight)法などがある。FS-BB法は、血液を低信号領域として黒く描出するBB (black blood)法において、ディフェーズ(diphase)傾斜磁場パルスとしてMPG (motion probing gradient)パルスを印加することによって関心領域(ROI: region of interest)における血液からの信号を選択的に低下させる手法である。このため、FS-BB法では、血液の位相が積極的に分散、すなわちdephasingされる。
従って、これらの高精度な位相補正が要求される条件下でイメージングを行う場合に、位相補正に誤差が存在するとアーチファクトが顕著に現れる。よって、撮像条件ごとにより特化した適切なAFI処理条件を設定することによって、アーチファクトの低減効果を向上させることができる。
次に、シミュレーションによって、適切なAFI処理条件のパラメータを決定する例を示す。まず、AFI処理の種類の決定例を示す。
図8は1次元のシミュレーションにより第1のAFI処理により生成したr-spaceデータI(x)の例を示す図、図9は1次元のシミュレーションにより第2のAFI処理により生成したr-spaceデータI(x)の例を示す図、図10は1次元のシミュレーションにより第3のAFI処理により生成したr-spaceデータI(x)の例を示す図、図11は1次元のシミュレーションにより非対称にサンプリングされたk-spaceに対して0-fillingを行って生成したr-spaceデータI0fill(x)の例を示す図、図12は図8から図12のシミュレーションに用いた部分的原k-spaceデータSorig(k)を示す図、図13は図12に示す部分的原k-spaceデータSorig(k)の作成に用いた全周波数領域におけるfull k-spaceデータSfull(k)を示す図、図14は図13に示すfull k-spaceデータSfull(k)を変換して得られるr-spaceデータI(x)及び位相分布Φ(x)を示す図である。
図8から図11において各横軸は、1次元のピクセル位置xを示し、各縦軸はピクセル位置xにおけるr-spaceデータI(x)を示す。また、図8から図11において、(A)はループ処理を行わない場合、つまりj=0の場合におけるr-spaceデータI0(x)を示し、(B)はループ処理を4回行った場合、つまりj=4の場合におけるr-spaceデータI4(x)を示す。
また、図12において横軸はk-spaceの位置kを示し、縦軸は位置kにおける部分的原k-spaceデータSorig(k)を示す。図12に示すようにKc=16とした部分的原k-spaceデータSorig(k) (-16≦k≦128)から図5、図6及び図7にそれぞれ示す第1、第2及び第3のAFI処理によって、それぞれ図8、図9及び図10に示すようなr-spaceデータI0(x), I4(x)を生成した。また比較参考用に図11の0-fillingを行って生成したr-spaceデータI0fill(x)を示す。
尚、図12に示す部分的原k-spaceデータSorig(k)は、図13に示す全周波数領域におけるfull k-spaceデータSfull(k)から非対称サンプリング領域-16≦k≦128のデータを切出して5位相分だけシフトさせたデータである。図14(A)は、図13に示す全周波数領域におけるfull k-spaceデータSfull(k)をFTして得られるr-spaceデータI(x)を示しており、図14(B)は、図14(A)に示すr-spaceデータI(x)から得られる位相分布Φ(x)を示している。
図8から図11に示すシミュレーション結果によれば、第2のAFI処理、すなわち低周波領域の位相分布Φlow(x)を用いた位相補正処理をhomodyne filter処理前に行う方法及び第3のAFI処理、すなわち全周波数領域の位相分布Φwhole(x)を用いた位相補正処理をhomodyne filter処理前に行う方法により、位相補正エラーが低減できることが分かる。
図5、図6及び図7に示す位相補正エラーを低減させるループ処理は、ループ処理前における位相補正のエラーが少ない程、収束が速くなる。従って、AFI処理の高速化の観点からは、第2のAFI処理又は第3のAFI処理が適していることが分かる。特に、部分的原k-spaceデータSorig(k)のピークがk-spaceの中心からシフトしている場合など位相誤差が大きい場合には、homodyne filter処理前に位相補正処理を行う第2のAFI処理又は第3のAFI処理を実行することが好適である。
また、ループ処理を行わない場合には、全周波数領域の位相分布Φwhole(x)を用いる第3のAFI処理が、最も位相補正エラーが少ないことが分かる。従って、データ処理時間を短くしつつ画質を確保する場合には、第3のAFI処理が適している。
また、第3のAFI処理において、ループ処理を実行すると、位相補正エラーが改善することが分かる。従って、一層の画質向上を重視する場合には、第3のAFI処理においてループ処理を実行することが適した処理条件となる。但し、シミュレーションによれば、ループ処理によって1次以下の低次の位相補正エラーが改善されるが、4回以上ループ処理を繰り返しても顕著な改善がないことが確認された。
一方、低周波領域の位相分布Φlow(x)を用いる第2のAFI処理では、ループ処理を実行しても位相補正エラーが殆ど改善されないことが分かる。従って、全周波数領域の位相分布Φwhole(x)が得られない場合には、データ処理時間を増加させないようループ処理を実行しない第2のAFI処理が適している。
逆に、位相補正処理をhomodyne filter処理後に行う第1のAFI処理の場合には、1次のピークシフト等の大きな位相変動があると位相補正によって十分に補正できないが、ループ処理の繰り返しによって位相補正エラーが改善されることが分かった。但し、第1のAFI処理の場合においても、4回以上ループ処理を繰り返しても顕著な改善がないことが確認された。
次に、シミュレーションによって、適切なhomodyne filterのゲインを決定する例を示す。
k-spaceデータの非対称サンプリング領域に対応するhomodyne filterの最大ゲインHhomo.maxの理論値は、位相補正後の位相補正エラーが無視できると仮定すると2≦Hhomo.max≦4となる。尚、従来は、非対称サンプリング領域に対応するhomodyne filterの最大ゲインHhomo.maxは固定値2であった。
しかし、位相分布の推定エラーが残存したり、推定エラーが位相補正によって強調される場合にhomodyne filterの最大ゲインHhomo.maxを2から4の値に設定すると、逆にアーチファクトが増加する場合がある。従って、homodyne filterの非対称サンプリング領域に対応する最大ゲインHhomo.maxは、理論値より小さい値が最適な場合がある。
homodyne filterの最大ゲインHhomo.max=1に設定することは、0-fillingを実行することに相当する。従って、homodyne filterの最大ゲインHhomo.max>1に設定すれば、少なくとも空間分解能が0-fillingより改善する。そこで、homodyne filterの最大ゲインを例えばHhomo.max=1.5に設定するなど、1<Hhomo.max<2に設定すれば、位相補正エラーを強調させずに空間分解能が0-fillingよりも改善されることが期待される。
homodyne filterの非対称サンプリング領域に対応する最大ゲインHhomo.maxは、第1、第2及び第3のAFI処理それぞれに適切な値に決定することができる。
homodyne filter処理より先に位相補正を行う第2及び第3のAFI処理では、位相補正後のr-spaceデータをk-spaceデータに変換すると非サンプリング部に信号が漏れる。すなわち、サンプリング部における信号成分が、非サンプリングにおける信号成分となる。
特に第3のAFI処理では全周波数領域の位相分布を用いた位相補正によって虚部信号がゼロに近い値となる。従ってk-spaceデータの実部はk=0に関して殆ど対称となる。このため、非サンプリング部分におけるk-spaceデータの信号強度は1/2程度となる。よって、第3のAFI処理を行う場合において位相エラーが無視できるのであれば、非サンプリング部分の信号強度がゼロになることによる画像のボケを改善するために、homodyne filterの最大ゲインHhomo.maxを2×2=4に設定することが妥当となる。
第2のAFI処理でも非サンプリング部に信号が漏れるため、位相エラーが無視できるのであれば、homodyne filterの最大ゲインHhomo.maxを2より大きな値に設定することが適切となる。
以上より、第1のAFI処理を行う場合、つまり位相補正処理をhomodyne filter処理後に行う場合には、位相補正エラーが強調されることを回避するために、homodyne filterの最大ゲインHhomo.maxを1<Hhomo.max≦2又は1<Hhomo.max<2となるように設定することが望ましい。また、第2及び第3のAFI処理を行う場合、つまり位相補正処理をhomodyne filter処理前に行う場合には、位相補正エラーを考慮しなければ、homodyne filterの最大ゲインHhomo.maxを2≦Hhomo.max≦4又は2<Hhomo.max≦4となるように設定することが望ましい。このように、homodyne filterの最大ゲインHhomo.maxを2以外の値に可変設定することができる。
次に、位相補正エラーの程度に着目してhomodyne filterの最大ゲインHhomo.maxを決定する例について説明する。
図15は、非対称サンプリング領域に対応するhomodyne filterのゲインHhomo.max=2とした第3のAFI処理のシミュレーション結果を示す図であり、図16は、非対称サンプリング領域に対応するhomodyne filterのゲインHhomo.max=4とした第3のAFI処理のシミュレーション結果を示す図である。
図15(A)は非対称サンプリング領域に対応する最大ゲインHhomo.max=2としたhomodyne filterのk方向における強度(ゲイン)Hhomoの分布を示す。位相補正後のk-spaceデータに対して図15(A)に示すhomodyne filterを適用し、ループ処理を伴わずに(j=0として)第3のAFI処理を実行すると図15(B)に示すr-spaceデータI0(x)が生成される。また、図15(A)に示すhomodyne filterを適用し、ループ処理を4回伴って(j=4として)第3のAFI処理を実行すると図15(C)に示すr-spaceデータI4(x)が生成される。
一方、図16(A)は非対称サンプリング領域に対応する最大ゲインHhomo.max=4としたhomodyne filterのk方向における強度(ゲイン)Hhomoの分布を示す。位相補正後のk-spaceデータに対して図16(A)に示すhomodyne filterを適用し、ループ処理を伴わずに(j=0として)第3のAFI処理を実行すると図16(B)に示すr-spaceデータI0(x)が生成される。また、図16(A)に示すhomodyne filterを適用し、ループ処理を4回伴って(j=4として)第3のAFI処理を実行すると図16(C)に示すr-spaceデータI4(x)が生成される。
図15及び図16によれば、ループ処理を行わない場合には、homodyne filterの最大ゲインHhomo.maxを4に設定すると、homodyne filterの最大ゲインHhomo.maxを2に設定した場合に比べて、高周波部分における信号強度が大きくなることが分かる。一方、ループ処理を行うと、homodyne filterの最大ゲインHhomo.maxを2に設定しても4に設定しても同様な信号強度分布のr-spaceデータI(x)が生成される。
これは、ループ処理において、r-spaceデータをk-spaceデータに変換した後、k-spaceデータの一部を原データと入れ換えるためであると考えられる。つまり、ループ処理によって位相補正後の信号の虚部がゼロとされ、位相補正のエラーが補正されるため、homodyne filterの最大ゲインHhomo.maxに依らず同等な信号強度分布のr-spaceデータI (x)が生成される。
また、homodyne filterの最大ゲインHhomo.max及びループ処理の回数jを変えて生成した各r-spaceデータI(x)のfull k-spaceデータSfull(k)に対する平均二乗誤差(RMSE: Root Mean Square Error)をそれぞれ計算した。その結果、Hhomo.max=2, j=0のときRMSE=1.44, Hhomo.max=2, j=4のときRMSE=1.55, Hhomo.max=4, j=0のときRMSE=1.98, Hhomo.max=4, j=4のときRMSE=1.54となった。すなわち、ループ処理を実行しない場合には、非対称サンプリング領域に対応するhomodyne filterの最大ゲインHhomo.max=2とした方が、RMSEを小さくすることができる。
従って、図15及び図16に示すシミュレーションに基づけば、空間分解能や位相補正エラーの条件が満たされていれば、RMSEを小さくする観点から非対称サンプリング領域に対応するhomodyne filterの最大ゲインHhomo.max=2に設定することが妥当となる。
図17は、第3のAFI処理において非対称サンプリング領域に対応するhomodyne filterの最大ゲインHhomo.maxを変えた場合におけるRMSEの変化を示す図であり、図18は、homodyne filterの強度分布形状を決定するパラメータを説明する図である。
図17において横軸は、homodyne filterの最大ゲインHhomo.maxを示し、縦軸は、-128≦x≦128の256 pixel 分のRMSEの総和を示す。また図18において横軸はk-spaceの位置kを示し、縦軸はhomodyne filterの強度Hhomoを示す。
第3のAFI処理、すなわち全周波数領域の位相分布Φwhole(x)を用いた位相補正後のk-space dataに、最大ゲインHhomo.maxを0.5から0.25刻みで4.0まで変化させてhomodyne filterを掛けた。そして、サンプリング領域の境界Kc=16とした部分的原k-spaceデータSorig(k) (-Kc=-16≦k≦Kmax=128)に対する第3のAFI処理によって生成された-128≦x≦128の全領域におけるr-space data I(x)のプロファイルとfull sampling data Sfull(k) (-Kmax=-128≦k≦Kmax=128)に対応する-128≦x≦128の全領域におけるr-spaceデータI(x)のプロファイルとの間においてRMSEを計算した。
図17の例では、homodyne filterの最大ゲインHhomo.max=1.75のときにRMSEが最小となる。RMSEが最小となるhomodyne filterの最大ゲインHhomo.maxの理論値は4である。従って、図17の例では、RMSEが最小となるhomodyne filterの最大ゲインHhomo.maxが理論値よりも大幅に小さい。この理由は、位相の空間分布に高周波成分が多いためであると考えられる。
以上より、第2及び第3のAFI処理では、非サンプリング部に信号が漏れることに着目すれば、homodyne filterの最大ゲインHhomo.maxを2≦Hhomo.max≦4とすることが妥当となる。
一方、位相の空間分布に高周波成分が多く位相補正エラーが無視できない場合には、RMSEを最小にするためにhomodyne filterの最大ゲインHhomo.maxを1<Hhomo.max<2とすることが妥当な場合がある。従って、位相補正に用いる位相の周波数分布に応じてhomodyne filterの最大ゲインHhomo.maxの取りうる範囲を1<Hhomo.max≦2又は2≦Hhomo.max≦4としてもよい。
尚、図17において、RMSEの計算に用いたhomodyne filterの強度分布形状Hhomoを決定するパラメータK1=16とした。すなわち、homodyne filterの強度分布形状Hhomoを決定するパラメータは、図18に示すように、サンプリング領域の境界Kcの他、K1及びK2がある。
サンプリング領域の境界Kcは、homodyne filterの強度Hhomoの立ち上がりの開始点となる。そして、境界Kcからk-space上の距離K1だけ進んだ位置においてhomodyne filterの強度Hhomoが一定の最大値Hhomo.maxとなる。また、homodyne filterの強度Hhomoが最大値Hhomo.maxの1/2となるときのk-space上の位置とhomodyne filterの強度Hhomoが最大値Hhomo.maxとなるときのk-space上の位置との間の距離はK2で表される。そして、K1の値が小さい程、フィルタの強度分布形状は矩形に近くなる。
図17の例では、r-space data I(x)のプロファイル全体の差に基づいてRMSEを計算したが、位相エラーの大きな領域のみのプロファイルの差に基づいてRMSEを計算してもよい。
非対称サンプリング領域に対応するhomodyne filterの最大ゲインHhomo.maxの最適値は、homodyne filterの対象となるデータに依存し、厳密には対称にサンプリングされたk-spaceデータがない限り計算することができない。そこで、上述したAFI処理の種類やRMSEに基づく方法の他、様々なアプローチで、homodyne filterの最大ゲインHhomo.maxを決定することができる。
例えば、サンプリング領域の境界Kcと同様にhomodyne filterの最大ゲインHhomo.maxを徐々に変化させてAFI処理によって複数の異なる最大ゲインHhomo.maxに対応する画像データを生成し、生成した複数の画像データを表示装置34に表示させるようにしてもよい。この場合、より画質が良好な画像データの選択情報を入力装置33からAFI処理条件設定部42に入力することによって、AFI処理条件設定部42は、画質が良好な画像データに対応するhomodyne filterの最大ゲインHhomo.maxをAFI処理条件として設定することが可能となる。画質が良好な画像データの選択は、ユーザの目視によって行うこともできるし、画像処理によって自動的に行うこともできる。
或いは、AFIに用いるシーケンスと所定のデータ収集条件が一致する同一種類のシーケンスを用いて同一の対象臓器又は同一の撮像部位から収集した他のMRデータに基づいて事前に実験的にhomodyne filterの最大ゲインHhomo.maxを決定することもできる。この場合、位相分布に影響を与えるTE等のデータ収集条件の一致度が、所要の割合であることが重要である。そして、事前に決定したhomodyne filterの最大ゲインHhomo.maxを以降の対応するAFIにおいて一定値として用いることができる。
更に別の方法としては、データ処理簡易化のために、AFI処理の種類に応じて取りうるhomodyne filterの最大ゲインHhomo.maxの範囲の平均値をAFI処理条件として決定することもできる。例えば、第1のAFI処理を行う場合には、homodyne filterの最大ゲインHhomo.maxを1<Hhomo.max≦2の中間値である1.5に固定し、第2又は第3のAFI処理を行う場合には、homodyne filterの最大ゲインHhomo.maxを2≦Hhomo.max≦4の中間値である3に固定することができる。
一方、適切なhomodyne filterの最大ゲインHhomo.maxを計算によって決定することもできる。ここでは、位相補正エラーの指標値を最小とする最適化計算によってhomodyne filterの最大ゲインHhomo.maxを計算する例と、位相補正によるk-spaceデータの非サンプリング部への信号の漏れの程度を表す指標値を用いた計算によってhomodyne filterの最大ゲインHhomo.maxを計算する例について説明する。
位相補正エラーの指標値は、例えば、部分的原k-spaceデータSorig(k)の0-fillingによって生成したr-spaceデータI0fill(x)とAFI処理によって生成したr-spaceデータIAFI(x)の差分値に用いて定義することができる。そして、homodyne filterの最大ゲインHhomo.maxが取り得る範囲内、例えば1<Hhomo.max≦2又は2≦Hhomo.max≦4において位相補正エラーの指標値を最小とする最適化計算によってhomodyne filterの最大ゲインHhomo.maxを決定することができる。
図19は、0-fillingによって生成したデータに基づく位相補正エラーの指標値が最小となるようにhomodyne filterの最大ゲインHhomo.maxを決定する方法を説明する図である。
図19において横軸はピクセル位置xを示し、縦軸はAFI処理によって生成されたr-spaceデータIAFI(x)と0-fillingによって生成したr-spaceデータI0fill(x)とのピクセル位置xにおける差の絶対値Isub(x)を示す。図19に示すように、AFI r-spaceデータIAFI(x)と0-filling r-spaceデータI0fill(x)の差の絶対値Isub(x)に対して閾値Ithを設定し、閾値Ithを超える絶対値Isub(x)の総和を位相補正エラーの指標値ErrorSumとすることができる。
すなわち、位相補正エラーの指標値ErrorSumは式(13)に示すアルゴリズムで計算することができる。
[数13]
for all pixels (x)
Isub(x) = |IAFI(x)-I0fill(x)|
if Isub(x) > Ith
ErrorSum = ErrorSum+Isub(x)
end
end
(13)
そして、位相補正エラーの指標値ErrorSumが最小となるように、homodyne filterの最大ゲインHhomo.maxを取りうる範囲内で変化させる収束計算を行うことによって、より適切なhomodyne filterの最大ゲインHhomo.maxを求めることができる。
また、0-filling r-spaceデータI0fill(x)とAFI r-spaceデータIAFI(x)との差に対する閾値処理により、差が大きい範囲を選択的に抽出して位相補正エラーの指標値ErrorSumの計算に用いることができる。これにより、0-filling r-spaceデータI0fill(x)とAFI r-spaceデータIAFI(x)との全ての位置における差を用いて位相補正エラーの指標値ErrorSumを計算する場合に生じるボケを回避することができる。
尚、式(13)で示すアルゴリズムでは、AFI r-spaceデータIAFI(x)と0-filling r-spaceデータI0fill(x)の振幅差の絶対値Isub(x)を加算した値を位相補正エラーの指標値ErrorSumとしたが、AFI r-spaceデータIAFI(x)と0-filling r-spaceデータI0fill(x)の位相差に基づいて位相補正エラーの指標値ErrorSumを定義してもよい。この場合にも位相差に対して閾値を設定することによりボケを低減させることができる。
他の例としては、0-filling r-spaceデータI0fill(x)の位相マップが所定の閾値を超える範囲に対応する領域をマスクしたAFI r-spaceデータIAFI(x)及び0-filling r-spaceデータI0fill(x)の2乗平均を位相補正エラーの指標値ErrorSumとして定義し、この指標値ErrorSumが最小となるような収束計算を行って最適なhomodyne filterの最大ゲインHhomo.maxを計算してもよい。
上述した位相補正エラーの指標値ErrorSumを最小にする収束計算は、画質向上のため、AFIを実行する度に毎回実行してもよい。但し、この収束計算は位相補正処理とhomodyne filter処理の繰り返しを伴うため、データ処理時間の増加に繋がる。そこで、AFIに用いるシーケンスと同一種類のシーケンスを用いて同一の対象臓器又は同一の撮像部位からデータを事前に収集し、同一種類のシーケンスで収集したデータに基づく収束計算によって位相補正エラーの指標値ErrorSumが最小となるhomodyne filterの最大ゲインHhomo.maxを求めてもよい。これにより、データ処理時間の増加を抑制することができる。
また、AFI又はAFIに先だって行うイメージングにおいてマルチスライスデータ収集を行う場合には、厳密にはスライスごとに収束計算の結果が異なる。そこで、複数のスライスに対応するhomodyne filterの最大ゲインHhomo.maxの平均値をAFI処理用のhomodyne filterの最大ゲインHhomo.maxとしたり、単一又は少数の代表スライスについてhomodyne filterの最大ゲインHhomo.maxを求めるようにすることもできる。これにより、データ処理時間の増加を抑制することができる。
尚、位相補正エラーを低減させるためのループ処理を行う場合には、ループ処理の繰り返しが、位相補正エラーの指標値ErrorSumが最小となるようにhomodyne filterの最大ゲインHhomo.maxを最適化することと等価である。従って、ループ処理を行う場合には、位相補正エラーの指標値ErrorSumが最小となるようにhomodyne filterの最大ゲインHhomo.maxの最適値を決定するプロセスは不要となる。
但し、ループ処理を行う場合であっても、homodyne filterの最大ゲインHhomo.maxが最適値に近い程、ループ処理の収束時間が短くなる。従って、ループ処理の時間短縮化のために、位相補正エラーの指標値ErrorSumを最小とするhomodyne filterの最大ゲインHhomo.maxの最適化を行ってもよい。
次に、第2及び第3のAFI処理において生じるk-spaceデータの非サンプリング部への信号の漏れの程度に対する指標を用いて適切なhomodyne filterの最大ゲインHhomo.maxを決定する方法について説明する。この方法によれば、位相補正エラーの指標値ErrorSumを用いる方法に比べて簡易に適切なhomodyne filterの最大ゲインHhomo.maxを決定することができる。
位相補正によるk-spaceにおける信号の漏れの程度に対する指標は、例えば、位相補正後におけるk-spaceデータSPcor(k)の信号の絶対値の対称度として定義することができる。
図20は、位相補正後におけるk-spaceデータSPcor(k)の信号の絶対値の例を示す図であり、図21は、図20に示す位相補正後のk-spaceデータSPcor(k)の生成に用いた元データを示す図である。
図20及び図21において各横軸はk-spaceの1D方向の位置kを示し、縦軸はk-spaceデータS(k)の信号値の絶対値を示す。
図20(A)は、第2のAFI処理における、つまり低周波領域の位相を用いた位相補正後のk-spaceデータSPcor_low(k)の信号の絶対値|SPcor_low(k)|のプロファイルを拡大した図であり、図20(B)は、第3のAFI処理における、つまり全周波領域の位相を用いた位相補正後のk-spaceデータSPcor_ whole(k)の信号の絶対値|SPcor_ whole(k)|のプロファイルを拡大した図である。
尚、図21(A)に示すfull sampling k-spaceデータSfull(k)の絶対値|Sfull(k)|からKc=16として-Kc≦k≦Kmaxの範囲を切り出した部分的原k-spaceデータSorig(k)の絶対値|Sorig(k)|が図21(B)である。そして、図21(B)に示す部分的原k-spaceデータSorig(k)の位相補正によって図20(A), (B)に示す位相補正後における各k-spaceデータSPcor_low(k), SPcor_ whole(k)が生成される。
図20(A)に示すように、低周波領域の位相を用いた位相補正を行うと、サンプリング領域の信号成分がk<-Kcの非サンプリング領域に若干漏れることが分かる。一方、図20(B)に示すように、全周波領域の位相を用いた位相補正を行うと、サンプリング領域の信号成分がk=0にほぼ対称となるようにk<-Kcの非サンプリング領域に漏れることが分かる。
図20(A), (B)に示すように、サンプリング領域から非サンプリング領域に漏れる信号成分の強度が大きくなる程、サンプリング領域における信号成分の強度が低下する。そこで、非サンプリング領域へ漏れる信号成分の強度が大きい程、homodyne filterの最大ゲインHhomo.maxが取り得る範囲内において大きくなるように、homodyne filterの最大ゲインHhomo.maxを設定すればよい。すなわち、非サンプリング領域に漏れる信号成分の強度に応じて適切なhomodyne filterの最大ゲインHhomo.maxを決定することができる。
ここでは、非サンプリング領域への信号の漏れの程度を表す指標として2つの定義例を示す。第1の指標SymIndex1は、位相補正前における部分的原k-spaceデータSorig(k)及び位相補正後におけるk-spaceデータSPcor(k)に基づいて定義することができる。すなわち、式(14)に示すように、第1の指標SymIndex1は、部分的原k-spaceデータSorig(k)のサンプリング領域の非対称部分Kc≦k≦Kmaxにおける信号強度の絶対値|Sorig(k)|の積算値の、位相補正後におけるk-spaceデータSPcor(k)のサンプリング領域の非対称部分Kc≦k≦Kmaxにおける信号強度の絶対値|SPcor(k)|の積算値に対する比として定義することができる。
式(14)のように第1の指標SymIndex1を定義すると、非サンプリング領域への信号の漏れが無ければサンプリング領域の非対称部分Kc≦k≦Kmaxにおける信号強度は位相補正前後において変化しないためSymIndex1=1となる。逆に、非サンプリング領域に信号がk=0に対称となるように漏れれば、位相補正後におけるk-spaceデータSPcor(k)のサンプリング領域の非対称部分Kc≦k≦Kmaxにおける信号強度は1/2になるためSymIndex1=2となる。
すなわち、1≦SymIndex1≦2である。従って、homodyne filterの最大ゲインHhomo.maxは、例えば式(15)により決定することができる。
[数15]
Hhomo.max=2a1/SymIndex1 (15)
但し、式(15)においてa1は0<a1≦1となる係数である。係数a1は、エラーの程度など条件により調整して決定することができる。
また、第2の指標SymIndex2は、位相補正後におけるk-spaceデータSPcor(k)のみに基づいて定義することができる。すなわち、式(16)に示すように、第2の指標SymIndex2は、位相補正後におけるk-spaceデータSPcor(k)のうち非サンプリング領域-Kmax≦k≦-Kcにおける信号強度の絶対値|SPcor(k)|の積算値の、サンプリング領域の非対称部分Kc≦k≦Kmaxにおける信号強度の絶対値|SPcor(k)|の積算値に対する比として定義することができる。
式(16)のように第2の指標SymIndex2を定義すると、非サンプリング領域への信号の漏れが無ければSymIndex2=0となる。逆に、非サンプリング領域に信号がk=0に対称となるように漏れれば、非サンプリング領域-Kmax≦k≦-Kcとサンプリング領域の非対称部分Kc≦k≦Kmaxにおいて位相補正後のk-spaceデータSPcor(k)は対称となるためSymIndex2=1となる。
すなわち、0≦SymIndex2≦1である。従って、homodyne filterの最大ゲインHhomo.maxは式(17)により決定することができる。
[数17]
Hhomo.max=2a2(1+SymIndex2) (17)
但し、式(17)においてa2は0<a2≦1となる係数である。係数a2は、エラーの程度など条件により調整して決定することができる。
ここまでは、homodyne filterの非対称サンプリング部分に対するゲインHhomo.maxを固定値として設定する場合について説明したが、一定でない値に設定することもできる。
図22はhomodyne filterの非対称サンプリング部分に対応するゲインを高周波部分ほど1に減衰するように設定した例を示す図である。
図22において横軸は、k-spaceにおける1D方向の位置kを示し、縦軸は位置kにおけるhomodyne filterのゲインHhomo(k)を示す。位相補正後におけるk-spaceデータSPcor(k)は、高周波領域ほど位相補正エラーが大きくなる可能性が高い。そこで、homodyne filterの非対称サンプリング部分に対応するゲインHhomo(k)に高周波部分ほど1に減衰する特性を与えることによって、位相補正エラーが大きくなる領域の信号強度を相対的に小さくすることができる。
尚、図22は、homodyne filterの最大ゲインHhomo.max(k)の取り得る範囲が1<Hhomo.max≦2である例を示しているが、homodyne filterの最大ゲインHhomo.max(k)の取り得る範囲が2≦Hhomo.max≦4である場合には、2≦Hhomo.max≦4の範囲で2に減衰するようにhomodyne filterのゲインHhomo(k)を設定すればよい。
以上のようにhomodyne filterの最大ゲインHhomo.max(k)を含む形状がAFI処理条件設定部42において設定される。この他、AFI処理に用いられる他のフィルタの強度(ゲイン)及び形状についてもAFI処理条件設定部42において設定される。
図23は、AFI処理において用いられるフィルタの具体的な形状の例を示す図であり、図24は、AFI処理において用いられるhomodyne filterの最大ゲインHhomo.max(k)をAFI処理の種類に応じて可変設定した例を示す図である。
図23(A), (B), (C)及び図24(A), (B)において、横軸は、k-spaceにおける1D方向の位置kを示す。
図23(A)は、非対称にサンプリングされた部分的原k-spaceデータSorig(k)から対称な低周波部分のk-spaceデータSlow(k) (-Kc≦k≦Kc)を抽出するためのLPF Hlow(k)の形状例を、図23(B)は、部分的原k-spaceデータSorig(k)から全周波数領域のデータ全体を切出すフィルタHwhole(k)の形状例を、図23(C)は、ループ処理内の合成処理のために部分的原k-spaceデータSorig(k)から-Kc≦k≦Kmaxの範囲を切り出すデータ切出しフィルタHmerge(k)の形状例を、それぞれ示す。
式(18-1)、式(18-2)、式(18-3)及び式(18-4)に示すように、それぞれLPF Hlow(k)、全周波数領域のデータ切出しフィルタHwhole(k)、合成処理用のデータ切出しフィルタHmerge(k)及びhomodyne filter Hhomo(k)を計算することができる。
[数18]
Hlow(k) = 1 :|k|≦Kc-K1
= exp[(-ln2){(k-(Kc-K1))/K2}2] : Kc-K1<|k|≦Kmax (18-1)
Hwhole(k) = Hlow(k) :k≦0
= 1 :0<k≦Kmax (18-2)
Hmerge(k) = Hlow(k) :k≦0
=1 :0<k≦Kmax (18-3)
Hhomo(k) = Hlow(k) :k≦0
= Hhomo.max-(Hhomo.max-1)Hlow(k) :0<k≦Kmax (18-4)
これにより、図23(A), (B)及び(C)に示すような形状のLPF Hlow(k)、全周波数領域のデータ切出しフィルタHwhole(k)及び合成処理用のデータ切出しフィルタHmerge(k)が生成される。尚、ゲインの遷移部分を滑らかに変化させるためにGauss関数が用いられている。
式(18-4)においてhomodyne filterの最大ゲインHhomo.max(k)は、AFI処理の種類に応じて可変設定することができる。例えば低周波領域の位相分布Φlow(x)を用いて位相補正を行う第1及び第2のAFI処理の場合には位相補正エラーを考慮してhomodyne filterの最大ゲインHhomo.max(k)=2とし、全周波数領域の位相分布Φwhole(x)を用いて位相補正を行う第3のAFI処理の場合には、homodyne filterの最大ゲインHhomo.max(k)=4とすることができる。
図24(A)は、homodyne filterの最大ゲインHhomo.max(k)=2とした第1及び第2のAFI処理用のhomodyne filterの形状Hhomo(k)を示す。また、図24(B)は、homodyne filterの最大ゲインHhomo.max(k)=4とした第3のAFI処理用のhomodyne filterの形状Hhomo(k)を示す。
式(18-1)、式(18-2)、式(18-3)及び式(18-4)並びに図23及び図24はk-spaceの負側の領域に非サンプリング領域を設定した場合の例である。k-spaceの正側の領域に非サンプリング領域を設定する場合には、k-spaceの負側と正側におけるフィルタ強度の値を反転すればよい。すなわち、図23及び図24において左右を反転すればよい。
尚、上述の説明では、AFI処理において、すなわちhomodyne filter Hhomo(k)又は全周波数領域のデータ切出しフィルタHwhole(k)によるフィルタ処理後のr-spaceデータIhomo(x)又はIwhole(x)に対して位相補正が行われる例を示したが、AFI処理の前処理としてk-spaceにおいて0次の位相補正及び1次の位相補正が行われるようにデータ処理条件を決定してもよい。
k-spaceデータS(k)に対する1次の位相補正はk-spaceデータのエコーピークの絶対値|Speak(k)|のk-spaceにおける位置を中心に合わせる位置シフトに相当する。従って、AFI処理の前処理として0次及び1次の位相補正又は少なくとも0次の位相補正を行えば、AFI処理における位相補正エラーの低減に繋げることができる。
特に、homodyne filter処理後に位相補正処理を行う第1のAFI処理の場合において、非対称の部分的原k-spaceデータSorig(k)k-spaceの位相誤差が大きいとAFI処理によって位相補正エラーが生じる。このため、第1のAFI処理を行う場合には、AFI処理の前処理として低次の位相補正を行うことが好適である。
また、式(19-1)に示すように、AFI処理によって生成されたAFI r-spaceデータIAFI(x)と0-fillingによって生成された0-filling r-spaceデータI0fill(x)とを重み係数wを用いて重み付け加算して得られるr-spaceデータIw(x)を診断画像データとしてもよい。この場合、重み係数wは、簡易な例として式(19-2)に示すようにAFI r-spaceデータIAFI(x)と0-filling r-spaceデータI0fill(x)との差とすることができる。
[数19]
Iw(x) = w*I0fill(x)+(1-w)* IAFI(x) (19-1)
w=|I0fill(x)- IAFI(x)| (19-2)
上述の0次及び1次の位相補正をAFI処理の前処理として行うか否か並びに式(19)に示す重み付け加算を行うか否かについてもAFI処理条件設定部42においてデータ処理条件として設定することができる。
次に磁気共鳴イメージング装置20の動作及び作用について説明する。
図25は、図1に示す磁気共鳴イメージング装置20によりAFI及びPIによって診断画像を収集する際の流れを示すフローチャートである。ここでは、非対称サンプリング領域の境界KCをプレスキャンによって決定する場合について述べる。
まずステップS1において、イメージングスキャンに先だってAFIにおけるデータサンプリング領域の境界Kcを決定するためのKc-PREPスキャンを含むプレスキャンが実行される。プレスキャンとしては、Kc-PREPスキャンの他、シミング用のプレスキャンやPI用にコイル要素24cの感度マップを測定するためのプレスキャンがある。
すなわち、撮像条件設定部40は、Kc-PREPスキャンを含むプレスキャン用の撮像条件を設定し、設定した撮像条件をシーケンスコントローラ31に出力する。これにより、後述するイメージングスキャンと同様な流れで必要なMRデータが収集される。シミング用のプレスキャン又はコイル要素24cの感度マップを測定するためのプレスキャンによって収集されたデータからは、AFI処理の位相補正に用いる位相分布を広範囲の周波数領域に亘って求めるともできる。
Kc-PREPスキャンにより収集されたMRデータからは、図4(B)に示すようにAFI処理によって画像データが生成される。従って、Kc-PREPスキャンに先だって、シミング用のプレスキャン等の他のプレスキャンを実行して位相分布を求めておけば、Kc-PREPスキャンにより収集されたデータに対する位相補正に利用することができる。そして、Kc-PREPスキャンにより収集された画像データから画質が良好な画像データを選択することにより、イメージングスキャン用のデータサンプリング領域の境界Kcを決定することができる。
次に、ステップS2において、k-spaceにおけるデータサンプリング領域-Kc≦k≦Kmaxに対応するMRデータの収集が実行される。すなわち、撮像条件設定部40は、非対称のデータサンプリング領域-Kc≦k≦KmaxからMRデータを収集するための所望のパルスシーケンスを含む撮像条件を設定する。撮像条件の設定は、表示装置34に表示された複数の像プロトコルの候補から入力装置33の操作によって所望の撮像プロトコルを選択することによって行うことができる。
そして、撮像条件設定部40は、パルスシーケンスを含む撮像条件をシーケンスコントローラ31に出力する。そうすると、シーケンスコントローラ31や静磁場用磁石21等のスキャンを実行するための磁気共鳴イメージング装置20の構成要素は、設定された撮像条件に従ってk-spaceの非対称サンプリング領域-Kc≦k≦Kmaxに対応するイメージングデータを被検体Pから収集する。
そのために、予め寝台37に被検体Pがセットされ、静磁場電源26により励磁された静磁場用磁石21(超伝導磁石)の撮像領域に静磁場が形成される。また、シミング用のプレスキャンによって収集されたデータに基づいて、シムコイル電源28からシムコイル22に所定の電流が供給されて撮像領域に形成された静磁場が均一化される。
そして、シーケンスコントローラ31は、撮像条件設定部40から受けたパルスシーケンスに従って傾斜磁場電源27、送信器29及び受信器30を駆動させることにより被検体Pがセットされた撮像領域に傾斜磁場を形成させるとともに、RFコイル24からRF信号を発生させる。
このため、被検体Pの内部における核磁気共鳴により生じたNMR信号が、RFコイル24により受信されて受信器30に与えられる。受信器30は、RFコイル24からNMR信号を受けて、デジタルデータのNMR信号である生データを生成する。受信器30は、生成した生データをシーケンスコントローラ31に与える。シーケンスコントローラ31は、生データをデータ処理部41に与え、データ処理部41は、k空間データベース43に形成されたk-spaceに生データを非対称サンプリング領域-Kc≦k≦Kmaxの部分的原k-spaceデータSorig(k)として配置する。
一方、ステップS3において、AFI処理条件設定部42は、撮像条件に対応する適切なAFI処理の種類、位相補正エラーを低減させるためのループ処理の回数j、homodyne filterの強度及び形状等のAFI処理のデータ処理条件を設定する。すなわち、位相補正処理をhomodyne filter処理前に行うかhomodyne filter処理後に行うか、全周波数領域の位相分布Φwhole(x)を位相補正に用いるか低周波領域の位相分布Φlow(x)を位相補正に用いるか、位相補正に用いる位相分布Φ(x)を補正対象となるk-spaceデータ自身のみから推定するか、他のデータのみから推定するか、自己データ及び他のデータの双方から推定するかといった事項が決定される。
例えば、AFI処理条件設定部42は、撮像条件設定部40から撮像プロトコルの選択情報を取得し、選択された撮像プロトコルに予め関連付けられたAFI処理条件をAFI条件データベース46から読み込むことによってAFI処理条件を設定することができる。また、homodyne filterの非対称サンプリング部に対応する最大ゲインHhomo.max(k)を部分的原k-spaceデータSorig(k)に基づく収束計算等の計算によって求める場合には、AFI処理条件設定部42において最大ゲインHhomo.max(k)が計算される。
また、AFI処理の種類が設定されると、AFI処理を構成する位相補正やフィルタ処理等の各種処理とPI処理との順序もAFI処理条件設定部42により決定される。例えば、処理時間短縮化のために、複数のコイル要素24cに対応する画像データの合成処理前に実行することが望ましいhomodyne filter処理等のAFI処理のうちの前段の処理の後にPI展開処理及び画像データの合成処理を実行し、合成処理後にAFI処理のループ処理等の残りの後段の処理が実行されるようにデータ処理順序が決定される。
更に、撮像条件に応じて必要であれば、AFI処理の前処理として0次及び1次の位相補正処理がデータ処理条件としてAFI処理条件設定部42において設定される。
そして、AFI処理条件設定部42は、前処理情報を含むAFI処理条件及びデータ処理順序をAFI処理部41Bに与える。
次に、ステップS4において、AFI処理部41Bは、AFI処理条件に従って、AFI処理の前処理及び前段の処理を実行する。処理後のデータは、AFI処理部41BからPI処理部41Aに与えられる。
次に、ステップS5において、PI処理部41Aは、PI展開処理及び画像データの合成処理を実行する。処理後のデータは、PI処理部41AからAFI処理部41Bに与えられる。
次に、ステップS6において、AFI処理部41Bは、AFI処理条件に従って、ループ処理等のAFI処理の後段の処理を実行する。これにより、対称サンプリングされたk-spaceデータの画像再構成処理によって生成された画像データと同等な画質の診断画像データが生成される。
次に、ステップS7において、データ処理部41は、PI処理及びAFI処理によって生成された診断用の画像データを表示装置34に表示させる。これにより、ユーザは、撮像条件に応じて適切な条件でAFI処理された高画質な診断画像を観察することが可能となる。
つまり以上のような磁気共鳴イメージング装置20は、AFIを行う際に、データ処理条件を撮像条件や位相分布に応じて適切に可変設定できるようにしたものである。
このため、磁気共鳴イメージング装置20によれば、AFIにおいて従来よりも画像の精度を向上させることができる。すなわち、k-spaceにおいて対称にサンプリングされたMRデータから生成された画像により近い画像を提供することができる。
また、位相補正処理をhomodyne filter処理前に行えば、エラーを低減させるためのループ処理の収束を速くできるため、画像精度向上のために必要となるデータ処理時間の増加を最小限に留めることができる。