JP7223619B2 - 磁気共鳴イメージング装置、画像処理装置及び画像処理方法 - Google Patents

Description

本発明の実施形態は磁気共鳴イメージング装置、画像処理装置及び画像処理方法に関する。

磁気共鳴イメージング装置は、静磁場中に置かれた被検体の原子核スピンをラーモア周波数の高周波（ＲＦ：Radio Frequency）信号で励起し、励起に伴って被検体から発生する磁気共鳴信号（ＭＲ（Magnetic Resonance）信号）を再構成して画像を生成する撮像装置である。

磁気共鳴イメージング技術の中に、Ｄｉｘｏｎ法と呼ばれる技術が知られている。Ｄｉｘｏｎ法では、異なるエコー時間ＴＥで収集したデータを再構成した複数の画像、例えば、第１のエコー時間に対応する第１画像と、第２のエコー時間に対応する第２画像とから、水と脂肪が分離された２つの画像、即ち、水画像（Ｗ画像）と脂肪画像（Ｆ画像）とを生成する技術である。水画像では、被検体の脂肪成分が抑制され、逆に、脂肪画像では、被検体の水成分が抑制されることになる。

第１画像及び第２画像の各画素値の位相は静磁場不均一性の影響を受ける。したがって、第１画像及び第２画像から水画像と脂肪画像とを算出する際には、静磁場不均一性の影響を正しく推定することが重要である。

そこで、静磁場不均一性の影響を正しく推定し、第１画像及び第２画像から水画像と脂肪画像とを正確に算出する技術が従来から種々提案されている。このような技術の中に、ＴＲＷ－Ｓ（sequential tree-reweighted message-passing：ツリー重み再配分メッセージ伝達法）アルゴリズムと呼ばれる技術がある。

ＴＲＷ－Ｓアルゴリズムは、静磁場が空間的に変化するとしても、その変化の仕方は連続的であり、かつ滑らかであろうという条件を利用して、静磁場不均一性の影響を正確に推定しようとするものである。

しかしながら、ＴＲＷ－Ｓアルゴリズムを用いても、なお、静磁場不均一性の影響を誤って推定する場合がある。静磁場不均一性の影響を誤って推定した場合、水画像の中の本来水成分であるべき領域の一部が誤って脂肪成分に入れ替わったり、逆に、脂肪画像の中の本来脂肪成分であるべき領域の一部が誤って水成分に入れ替わったりするという事象が起こり得る。

Johan Berglund et al., Two-point Dixon Method With Flexible Echo Times, MRM 65:994-1004 (2011)

本発明が解決しようとする課題は、ＴＲＷ－Ｓアルゴリズムを用いた水画像と脂肪画像の分離技術において、水成分の領域と脂肪成分の領域とを、誤りが少なく、より正確に推定できるようにすることである。

実施形態の磁気共鳴イメージング装置は、第１のエコー時間で収集されるデータから第１画像を再構成し、第２のエコー時間で収集されるデータから第２画像を再構成する再構成部と、前記第１画像と前記第２画像とから、画素ごとに、一方が真であり他方が偽である第１位相及び第２位相を算出し、第１位相マップ及び第２位相マップを夫々生成する位相算出部と、前記第１位相マップ及び第２位相マップから、ＴＲＷ－Ｓ（sequential tree-reweighted message-passing）アルゴリズムを用いて、被検体の撮像時の静磁場分布に対応する真の位相分布を推定した推定真位相マップを生成するＴＲＷ－Ｓ処理部と、
前記第１及び第２画像と推定真位相マップとを用いて、水画像と脂肪画像とを生成する画像生成部と、を備え、前記ＴＲＷ－Ｓ処理部は、（ａ）前記推定真位相マップの誤推定領域を検出し、（ｂ）検出した前記誤推定領域が前記第１位相であるべき領域を前記第２位相の領域であると誤って推定している場合には、第２位相マップの該当領域を前記第１位相に置換し、検出した前記誤推定領域が前記第２位相であるべき領域を前記第１位相の領域であると誤って推定している場合には、前記第１位相マップの該当領域を前記第２位相に置換し、（ｃ）置換した領域の位相を固定した状態で、前記ＴＲＷ－Ｓアルゴリズムを実行して、前記推定真位相マップを更新し、前記画像生成部は、更新された前記推定真位相マップと、前記第１及び第２画像とを用いて、前記水画像と前記脂肪画像とを生成する。

第１の実施形態の磁気共鳴イメージング装置の構成例を示す図。 本実施形態の磁気共鳴イメージング装置で行っているＤｉｘｏｎ法の全体的な概要を示す図。 水領域と脂肪領域の入れ替わりが起こったときの水画像と脂肪画像を模式的に例示する図。 改善型ＴＲＷ－Ｓアルゴリズムを実行する磁気共鳴イメージング装置の機能構成例を示すブロック図。 改善型ＴＲＷ－Ｓアルゴリズムを実行する磁気共鳴イメージング装置の動作例を示す全体フローチャート。 従来のＴＲＷ－Ｓアルゴリズムの課題を説明する図。 改善型ＴＲＷ－Ｓ機能の詳細動作例を示す第１の個別フローチャート。 改善型ＴＲＷ－Ｓ機能の詳細動作例を示す第２の個別フローチャート。 本実施形態の改善型ＴＲＷ－Ｓアルゴリズムの第１の説明図。 本実施形態の改善型ＴＲＷ－Ｓアルゴリズムの第２の説明図。 本実施形態の改善型ＴＲＷ－Ｓアルゴリズムの第３の説明図。 本実施形態の改善型ＴＲＷ－Ｓアルゴリズムの第４の説明図。 改善型ＴＲＷ－Ｓアルゴリズムを実行する画像処理装置の機能構成例を示すブロック図。

実施形態に係る磁気共鳴イメージング装置及び画像処理装置について、添付図面を参照して説明する。なお、以下の実施形態では、同一の参照符号を付した部分は同様の動作をするものとして、重複する説明を適宜省略する。
（磁気共鳴イメージング装置）

図１は、実施形態に係る磁気共鳴イメージング装置１の全体構成を示すブロック図である。磁気共鳴イメージング装置１は、磁石架台１００、制御キャビネット３００、コンソール４００、寝台５００、及びＲＦ（Radio Frequency）コイル２０を備える。

磁石架台１００は、静磁場磁石１０、傾斜磁場コイル１１、及びＷＢ（Whole Body）コイル１２を有している。これらの構成品は円筒状の筐体に収納されている。寝台５００は、寝台本体５０と天板５１を有している。

制御キャビネット３００は、傾斜磁場電源３１（Ｘ軸用３１ｘ、Ｙ軸用３１ｙ、Ｚ軸用３１ｚ）、コイル選択回路３６、ＲＦ受信器３２、ＲＦ送信器３３、及びシーケンスコントローラ３４を備えている。

コンソール４００は、処理回路４０、記憶回路４１、ディスプレイ４２、及び入力デバイス４３を備えている。コンソール４００は、ホスト計算機として機能する。

磁石架台１００の静磁場磁石１０は、概略円筒形状をなしており、被検体、例えば患者、が搬送されるボア内に静磁場を発生させる。ボアとは、磁石架台１００の円筒内部の空間のことである。静磁場磁石１０は超電導コイルを内蔵し、液体ヘリウムによって超電導コイルが極低温に冷却されている。静磁場磁石１０は、励磁モードにおいて静磁場用電源（図示せず）から供給される電流を超電導コイルに印加することで静磁場を発生する。その後、永久電流モードに移行すると、静磁場用電源は切り離される。一旦永久電流モードに移行すると、静磁場磁石１０は長時間、例えば１年以上に亘って、大きな静磁場を発生し続ける。

傾斜磁場コイル１１も概略円筒形状をなし、静磁場磁石１０の内側に固定されている。この傾斜磁場コイル１１は、傾斜磁場電源（３１ｘ、３１ｙ、３１ｚ）から供給される電流によりＸ軸，Ｙ軸，Ｚ軸の方向に傾斜磁場を被検体に印加する。

寝台５００の寝台本体５０は天板５１を上下方向及び水平方向に移動することができる。撮像前に天板５１に載置された被検体を所定の高さまで移動させる。その後、撮像時には天板５１を水平方向に移動させて被検体をボア内に移動させる。

ＷＢコイル１２は全身用コイルとも呼ばれ、傾斜磁場コイル１１の内側に被検体を取り囲むように概略円筒形状に固定されている。ＷＢコイル１２は、ＲＦ送信器３３から伝送されるＲＦパルスを被検体に向けて送信する。また、水素原子核の励起によって被検体から放出される磁気共鳴信号、即ちＭＲ（Magnetic Resonance）信号を受信する。

磁気共鳴イメージング装置１は、ＷＢコイル１２の他、図１に示すようにＲＦコイル２０を備える。ＲＦコイル２０は、被検体の体表面に近接して載置されるコイルである。ＲＦコイル２０は複数の要素コイルを備えている。これら複数の要素コイルは、ＲＦコイル２０の内部でアレイ状に配列されるため、ＰＡＣ（Phased Array Coil）と呼ばれることもある。ＲＦコイル２０には幾つかの種別がある。例えば、ＲＦコイル２０として、図１に示すように被検体の胸部や腹部、或いは脚部に設置されるボディコイル（Body Coil）や、被検体の背側に設置されるスパインコイル（Spine Coil）といった種別がある。

ＲＦ送信器３３は、シーケンスコントローラ３４からの指示に基づいてＲＦパルスを生成する。生成したＲＦパルスはＷＢコイル１２に伝送され、被検体に印加される。ＲＦパルスの印加によって被検体からＭＲ信号が発生する。このＭＲ信号をＲＦコイル２０又はＷＢコイル１１が受信する。

ＲＦコイル２０で受信したＭＲ信号、より具体的には、ＲＦコイル２０内の各要素コイルで受信したＭＲ信号は、天板５１及び寝台本体５０に設けられたケーブルを介してコイル選択回路３６に伝送される。コイル選択回路３６は、ＲＦコイル２０から出力される信号、或いはＷＢコイルから出力される信号を、シーケンスコントローラ３４或いはコンソール４００から出力される制御信号に応じて選択する。

選択された信号はＲＦ受信器３２に出力される。ＲＦ受信器３２は、チャネル信号、即ちＭＲ信号をＡＤ（Analog to Digital）変換して、シーケンスコントローラ３４に出力する。デジタルに変化されたＭＲ信号は、生データ（Raw Data）と呼ばれることもある。なお、ＡＤ変換は、ＲＦコイル２０の内部やコイル選択回路３６で行ってもよい。

シーケンスコントローラ３４は、コンソール４００による制御のもと、傾斜磁場電源３１、ＲＦ送信器３３及びＲＦ受信器３２をそれぞれ駆動することによって被検体のスキャンを行う。スキャンによってＲＦ受信器３２から生データを受信すると、シーケンスコントローラ３４は、その生データをコンソール４００に送信する。

シーケンスコントローラ３４は、処理回路（図示を省略）を具備している。この処理回路は、例えば所定のプログラムを実行するプロセッサや、ＦＰＧＡ（Field Programmable Gate Array）、ＡＳＩＣ（Application Specific Integrated Circuit）等のハードウェアで構成される。

コンソール４００は、記憶回路４１、入力デバイス４３、ディスプレイ４２、及び処理回路４０を備える。記憶回路４１は、ＲＯＭ（Read Only Memory）やＲＡＭ（Random Access Memory）の他、ＨＤＤ（Hard Disk Drive）や光ディスク装置等の外部記憶装置を含む記憶媒体である。記憶回路４１は、各種の情報やデータを記憶する他、処理回路４０が具備するプロセッサが実行する各種のプログラムを記憶する。

入力デバイス４３は、例えば、マウス、キーボード、トラックボール、タッチパネル等であり、各種の情報やデータを操作者が入力するための種々のデバイスを含む。ディスプレイ４２は、液晶ディスプレイパネル、プラズマディスプレイパネル、有機ＥＬパネル等の表示デバイスである。

処理回路４０は、例えば、ＣＰＵや、専用又は汎用のプロセッサを備える回路である。プロセッサは、記憶回路４１に記憶した各種のプログラムを実行することによって、後述する各種の機能を実現する。処理回路４０は、ＦＰＧＡやＡＳＩＣ等のハードウェアで構成してもよい。これらのハードウェアによっても後述する各種の機能を実現することができる。また、処理回路４０は、プロセッサとプログラムによるソフトウェア処理と、ハードウェア処理とを組わせて、各種の機能を実現することもできる。

なお、図１に示す磁気共鳴イメージング装置１の構成のうち、コンソール４００以外の構成品（制御キャビネット３００、磁石架台１００及び寝台５００）で、撮像部を構成している。撮像部は、スキャナと呼ぶ場合がある。

実施形態の磁気共鳴イメージング装置１では、Ｄｉｘｏｎ法と呼ばれる撮像法、或いは、画像処理方法を用いている。図２は、本実施形態の磁気共鳴イメージング装置１で行っているＤｉｘｏｎ法の全体的な概要を示す図である。

Ｄｉｘｏｎ法では、異なるエコー時間ＴＥで収集したデータを再構成した複数の画像、例えば、第１のエコー時間（ＴＥ１）に対応する第１画像と、第２のエコー時間（ＴＥ２）に対応する第２画像とから、水と脂肪が分離された２つの画像、即ち、水画像（Ｗ画像）と脂肪画像（Ｆ画像）とを生成する技術である。水画像では、被検体の脂肪成分が抑制され、逆に、脂肪画像では、被検体の水成分が抑制されることになる。なお、以下では、第1画像と第２画像の２つの画像から水画像と脂肪画像を生成する、所謂、２-ポイントＤｉｘｏｎ法に基づいて説明するが、本実施形態はこれに限定されるものではなく、３以上の複数のエコー時間に対応する複数の画像を用いて水画像と脂肪画像を生成する、所謂、マルチポイントＤｉｘｏｎ法に拡張することもできる。

第１画像及び第２画像の各画素値の位相は静磁場不均一性の影響を受ける。したがって、第１画像及び第２画像から水画像と脂肪画像とを算出する際には、静磁場不均一性の影響を正しく推定すること重要である。

第１画像の複素画素値Ｓ_１と第２画像の複素画素値Ｓ_２とは、夫々、静磁場不均一性ΔＢに起因する位相φの影響を受けている。後述するように、位相φは、複素画素値Ｓ_１及び第２画像の複素画素値Ｓ_２、並びに、夫々の画像に対応するエコー時間ＴＥ１、ＴＥ２等の既知のパラメータを含む方程式を解くことによって求めることができる。しかしながら、この方程式は２次方程式であるため、この２次方程式を解くと、正と負の符号に夫々対応する２つの位相の解、φ_Ａとφ_Ｂとが求まることになる。

位相φ_Ａと位相φ_Ｂは、互いにペアとして算出され、一方が真の位相であり、他方が偽の位相を表わしている。位相φ_Ａと位相φ_Ｂのペアは、第１、第２画像の画素毎に算出されるため、位相φ_Ａと位相φ_Ｂの集合は、これらを対応する画素位置に配列することにより、図２の中央上段に示すように、位相マップφ_Ａと位相マップφ_Ｂとして表すことができる。

静磁場不均一性を正しく判定するためには、位相φ_Ａと位相φ_Ｂのペアのうち、どちらが真でどちらが偽であるかを正しく判定し、位相φ_Ａと位相φ_Ｂのどちらか一方を真の位相φとして推定する必要がある。

ＴＲＷ－Ｓアルゴリズムは、静磁場が空間的に変化するとしても、その変化の仕方は連続的であり、かつ滑らかであろうという条件を利用して、静磁場不均一性の影響を正確に推定しようとするものである。非特許文献１は、従来のＴＲＷ－Ｓアルゴリズムについて詳細に開示している。非特許文献１では、ＴＲＷ－Ｓアルゴリズムの実施に先立って、処理負荷の低減等のために、位相マップφ_Ａと位相マップφ_Ｂの解像度を低解像度化する処理を行っている。そして、図２の右部に示すように、低解像度化した位相マップΦ_Ａと位相マップΦ_Ｂに対してＴＲＷ－Ｓアルゴリズムを適用し、真の位相マップを得るものとしている。

しかしながら、非特許文献１に記載のＴＲＷ－Ｓアルゴリズムを適用したとしても、なお、真の位相を誤って推定する場合がある。即ち、位相Φ_Ａと位相Φ_Ｂのうち、誤って偽の方の位相を選択する場合がある。そして、偽の位相を選択し、その位相を用いて水画像と脂肪画像とを生成すると、水画像の中の本来水成分であるべき領域の一部が誤って脂肪成分に入れ替わったり、逆に、脂肪画像の中の本来脂肪成分であるべき領域の一部が誤って水成分に入れ替わったりするという事象が起こり得る。

図３は、水領域と脂肪領域の入れ替わりが起こったときの、水画像と脂肪画像を模式的に例示する図である。図３は、肝臓を含むアキシャル断面の一例を示しており、図３の上段は、水領域と脂肪領域の入れ替わりが発生していない、正しい水画像と脂肪画像とを例示している。一方、図３の下段は、水領域と脂肪領域の入れ替わりが発生し、その結果、誤って生成された水画像と脂肪画像とを例示している。

本実施形態の磁気共鳴イメージング装置１は、このような水領域と脂肪領域の入れ替わり事象に対処するものであり、特に、従来のＴＲＷ－Ｓアルゴリズムを改善することにより、水領域と脂肪領域の入れ替わり事象を抑制するものである。

図４は、改善型ＴＲＷ－Ｓアルゴリズムを実行するＤｉｘｏｎ法の処理機能に焦点を絞った処理回路４０の機能構成例を示すブロック図である。また、図５は、改善型ＴＲＷ－Ｓアルゴリズムを実行するＤｉｘｏｎ法の処理機能に焦点を絞った磁気共鳴イメージング装置１の動作例を示す全体フローチャートである。

図４に示すように、コンソール４００の処理回路４０は、撮像条件設定機能４０１、再構成機能４０２、第１画像（ＴＥ＝ＴＥ１）生成機能４０３、第２画像（ＴＥ＝ＴＥ２）生成機能４０４、第１位相算出機能４０５、第２位相算出機能４０６、改善型ＴＲＷ－Ｓ機能４０７、及び、画像生成機能４０８を実現する。これらの各機能は、例えば、処理回路４０が具備するプロセッサが、記憶回路４１に記憶される所定のプログラムを実行することによって実現される。これら各機能の動作について図５の全体フローチャートを用いて説明する。また、特に、改善型ＴＲＷ－Ｓ機能４０７に関しては、図７及図８に示す個別フローチャートと、関連する幾つかの図面を用いて説明する。

まず、図５のステップＳＴ１０において第１画像を生成する。また、ステップＳＴ２０において第２画像を生成する。具体的には、図４に示す撮像条件設定機能４０１が、第１画像用の第１のパルスシーケンスを設定し、この第１のパルスシーケンスの実行によって収集されたＭＲ信号を、第１画像生成機能４０３が再構成することにより、第１画像を生成する。同様に、図４に示す撮像条件設定機能４０１が、第２画像用の第２のパルスシーケンスを設定し、この第２のパルスシーケンスの実行によって収集されたＭＲ信号を、第２画像生成機能４０４が再構成することにより、第２画像を生成する。

第１と第２のパルスシーケンスは、エコー時間ＴＥが異なる値に設定されるものの、パルスシーケンスの種類や、エコー時間以外のパルスシーケンスのパラメータは、基本的に同一に設定される。第１、第２のパルスシーケンスは、ＳＥ（Spin Echo）系のパルスシーケンスでもよいし、ＧＲＥ（Gradient Echo）系のパルスシーケンスでもよい。

一方、第１のパルスシーケンスのエコー時間はＴＥ１に設定され、第２のパルスシーケンスのエコー時間はＴＥ２に設定される。周知のように、水成分と脂肪成分とでは、磁気共鳴周波数が異なっているため（例えば、互いに３．５ｐｐｍ程度離れているため）、エコー時間ＴＥが異なると、水成分の位相の回転量と、脂肪成分の位相の回転量とが異なる。水成分と脂肪成分の位相が逆相となるようなＴＥ１とＴＥ２の設定によって収集されたＭＲ信号から生成された画像はOut-of-Phase画像と呼ばれる。また、水成分と脂肪成分の位相が同相となるようなＴＥ１とＴＥ２の設定によって収集されたＭＲ信号から生成された画像はIn-Phase画像と呼ばれる。

水画像と脂肪画像を生成するためには、エコー時間ＴＥ１とＴＥ２が既知であれば十分であり、第１画像と第２画像の一方がOut-of-Phase画像であり、他方がIn-Phase画像である必要はない。ただし、以下では、便宜上、第１画像がOut-of-Phase画像であり、第２画像がIn-Phase画像であるものとして説明する。なお、以下の説明で用いている数式は、主に非特許文献１に記載されている数式に基づいている。

第１画像と第２画像とから、位相φ_Ａと位相φ_Ｂを算出する過程を以下に説明する。第１画像の画素値（複素数）をＳ_１とし、第２画像の画素値（複素数）をＳ_２とすると、画素値Ｓ_１とＳ_２を以下の式で表すことができる。

（式１）、（式２）において、Ｗ及びＦは、夫々、水と脂肪のエコー時間ＴＥ＝０の時の信号値（実数）である。

ｂ_０は、ＴＥ＝ＴＥ１のときの水信号の位相をφ_０とするときの位相項（ｂ_０＝exp（ｊφ_０））である。また、ｂは、ＴＥ２とＴＥ１との間の時間ΔＴ（＝ＴＥ２－ＴＥ１）の間に、主に静磁場不均一性に起因して回転する位相をφとするときの位相項（ｂ＝exp（ｊφ））である。パラメータａ_１とａ_２は、以下の（式３）、（式４）で表される。

上式において、γは磁気回転比、Ｂ_０は静磁場強度、δは水に対する脂肪のケミカルシフト（約３．５ｐｐｍ）であり、ＴＥ１とＴＥ２は、励起パルス印加時をＴＥ＝０とするときの、第１画像、第２画像の収集時のエコー時間である。このように、パラメータａ_１とａ_２は、撮像条件として設定したエコー時間ＴＥ等により定まる既知のパラメータである。

（式１）、（式２）を夫々２乗すると、以下の（式５）、（式６）が得られる。

ここで、脂肪比率として、Ｑ＝Ｆ／（Ｗ＋Ｆ）を導入すると、（式５）、（式６）は、Ｑの２次方程式となる。そして、この２次方程式の解は、次の（式７）となる。

（式７）において、ｃ１、ｃ２、ｃ３は、以下の（式８）－（式１０）である。

が得られる。ここで、（式１１）、（式１２）のｂ_Ａ、ｂ_Ｂは、位相φ_Ａと位相φ_Ｂと以下の（式１３）、（式１４）で関係づけられている。
ｂ_Ａ＝ｅｘｐ（ｊφ_Ａ） （式１３）
ｂ_Ｂ＝ｅｘｐ（ｊφ_Ｂ） （式１４）

上記のように、（式１）から（式１４）に基づいて、第１、第２画像の画素値Ｓ_１、Ｓ_２、及び、既知のパラメータａ_１、ａ_２から、位相φ_Ａと位相φ_Ｂを算出することができる。ｂ_Ａとｂ_Ｂは（式７）における±の符号の夫々に対応するものであり、同様に、位相φ_Ａと位相φ_Ｂも（式７）における±の符号の夫々に対応するものである。位相φ_Ａと位相φ_Ｂのうち、一方が実際の静磁場分布に対応する真の位相であり、他方が偽の位相となる。

（式１）から（式１４）までの演算処理が、図５のステップＳＴ３０の処理に該当する。位相φ_Ａと位相φ_Ｂは第１、第２画像の画素値Ｓ_１、Ｓ_２毎に算出される。そして、位相φ_Ａと位相φ_Ｂを該当する画素位置に配列したものが、ステップＳＴ４０で生成される位相マップφ_Ａと位相マップφ_Ｂである。

ステップＳＴ５０では、第１、第２画像の画素値Ｓ_１、Ｓ_２から、例えば、以下の式（１５）により、振幅ｍを算出し、これを該当する画素位置に配列した振幅マップを生成する。
ｍ＝｜Ｓ_１｜＋｜Ｓ_２｜ （式１５）
なお、（式１５）に換えて、In-Phase画像（第２画像）の画素値の振幅のみを用いて振幅マップを生成してもよい。

ステップＳＴ６０では、ステップＳＴ４０で生成した位相マップφ_Ａと位相マップφ_Ｂを低解像化して、低解像位相マップΦ_Ａと位相マップΦ_Ｂを夫々生成する。ステップＳＴ８０で行う改善型ＴＲＷ－Ｓアルゴリズムの処理負荷を低減するためである。同様に、ステップＳＴ７０では、ステップＳＴ５０で生成した振幅マップ（ｍ）を低解像化して、低解像振幅マップ（Ｍ）を生成する。

ステップＳＴ８０で、低解像位相マップΦ_Ａ、Φ_Ｂ、及び、低解像振幅マップ（Ｍ）を用いて改善型ＴＲＷ－Ｓを実行し、低解像位相マップΦ_Ａ、Φ_Ｂから真の位相Φを推定して、低解像度の推定真位相マップ（Φ）を生成する。ステップＳＴ８０の具体的な処理に関しては、後段で詳述する。

ステップＳＴ９０では、ステップＳＴ８０で生成した低解像度の推定真位相マップ（Φ）をスムージングして、原画像（第１、第２画像）の解像度に高解像度化した推定真位相マップ（φ）を生成する。

最後に、ステップＳＴ１００では、第１、第２画像、及び、推定真位相マップ（φ）を用いて、以下の各式に基づいて、水画像と脂肪画像の夫々の画素値ＷとＦを生成する。

（式１６）内の「ｂ」が、改善型ＴＲＷ－Ｓアルゴリズムによって推定される真の位相φに対応する位相項（ｂ＝ｅｘｐ（ｊφ））である。（式１９）のＷ_ＬＳ、Ｆ_ＬＳが、最終的に得られる水画像と脂肪画像の画素値である。なお、（式２０）は、Ｗ_ＬＳ、Ｆ_ＬＳの二乗誤差を最少とするための逆行列（ムーア・ペンローズ逆行列）である。

（改善型ＴＲＷ－Ｓアルゴリズム）
以下、ステップＳＴ８０の改善型ＴＲＷ－Ｓアルゴリズムについて説明する。なお、本実施形態の改善型ＴＲＷ－Ｓアルゴリズムは解像度の高低に依存しない。したがって、以下では、ステップＳＴ４０とステップＳＴ５０で生成した、低解像化処理前の位相マップφ_Ａと位相マップφ_Ｂ、及び、振幅マップ（ｍ）を用いて、本実施形態の改善型ＴＲＷ－Ｓアルゴリズムを説明する。

改善型ＴＲＷ－Ｓアルゴリズムの説明に先立ち、従来のＴＲＷ－Ｓアルゴリズムの課題について、図６を用いて説明する。図６の左側は、ステップＳＴ４０で生成した位相マップφ_Ａと位相マップφ_Ｂを例示する図である。以下の説明では、位相を大きさ１のベクトルで表記するものとしている。さらに、白色の矢印で表記したベクトルによって位相φ_Ａを表わすものとし、黒色の矢印で表記したベクトルによって位相φ_Ｂを表わすものとしている。

図６の右側が、位相マップφ_Ａ及び位相マップφ_Ｂから、ＴＲＷ－Ｓアルゴリズムによって推定された真の位相マップ（φ）（以下、推定真位相マップ（φ）と呼ぶ）である。位相φは、静磁場分布に対応するものであるため、位相φの空間分布は本来連続的な分布を示すはずである。しかしながら、ＴＲＷ－Ｓアルゴリズムによって得られた推定真位相マップ（φ）では、左上の領域のベクトル（黒いベクトル）が、それ以外の領域のベクトル（白いベクトル）とほぼ逆向きになっており、両者の領域の境界で不連続となっている。したがって、左上の黒いベクトルの領域は、位相φを誤って推定した可能性が高い。

改善型ＴＲＷ－Ｓアルゴリズムでは、このような誤推定を可能な限り低減することを狙っている。図７及び図８は、改善型ＴＲＷ－Ｓアルゴリズム（図５のステップＳＴ８０）の動作例を示すフローチャートである。以下、このフローチャートと、図９乃至図１２を用いて、本実施形態の改善型ＴＲＷ－Ｓアルゴリズムを説明する。

ステップＳＴ２００では、従来のベーシックＴＲＷ－Ｓ処理（非特許文献１に記載のＴＲＷ－Ｓ処理）を行う。ベーシックＴＲＷ－Ｓ処理では、位相マップφ_Ａと位相マップφ_Ｂを取得し、推定真位相マップ（φ）を生成する。ベーシックＴＲＷ－Ｓ処理では、推定真位相マップ（φ）を繰り返し演算によって求めている。そして、最終的な推定真位相マップ（φ）の生成の過程において、途中段階の推定真位相マップ（φ）に対して、以下の（式２１）で表される値をコスト画素とするマップ（以下、コストマップと呼ぶ）を生成している。

（式２１）において、ｂｓ、ｂｔは、推定真位相マップ（φ）内において、隣接する位相φ_ｓとφ_ｔに対応する位相項である。また、ｍ_ｓ、ｍ_ｔは、振幅マップ（ｍ）内において隣接する振幅値である。

また、コストマップの各コスト画素を積算して、以下の（式２２）で表れるコスト値Ｃを算出している。

ベーシックＴＲＷ－Ｓでは、（式２２）で示されるコスト値が所定値以下となった時点で処理を終了し、その時の推定真位相マップを、最終的な推定真位相マップとしている。

これに対して、本実施形態では、ステップＳＴ２００のベーシックＴＲＷ－Ｓに引き続き、以下に説明するステップＳＴ２０１からステップＳＴ２０６の処理を行っている。

ステップＳＴ２０１では、処理方向と、処理方向の順逆の初期設定を行う。そして、ステップＳＴ２０２からステップＳＴ２０４のループにより、処理方向と処理方向の順逆を更新しながら、ステップＳＴ２０３の処理を繰り返し行う。

ここで、処理方向の設定とは、互いに直交する３方向（例えば、物理座標系におけるＸ方向、Ｙ方向、Ｚ方向や、３次元画像データにおける、行方向、列方向、厚み方向等）のうち、どの方向の処理を行うかを設定することである。また、処理方向の順逆の設定とは、設定された処理方向において、順方向（正の方向）と逆方向（負の方向）のどちらの向きの処理を行うかを設定することである。したがって、第１、第２画像が３次元画像の場合、直交３方向と夫々の方向の順方向と逆方向の合計６つの方向に対して、ステップＳＴ２０３の変形ＴＲＷ－Ｓ処理を実施することになる。

図８は、ステップＳＴ２０３のより詳細な処理例を示すフローチャートである。図８のステップＳＴ３００では、推定真位相マップ内の不連続位置をコスト壁として検出する処理を行う。図９は、コスト壁検出の概念を説明する図である。

図９（ａ）は、図６の右図と同じ図であり、ベーシックＴＲＷ－Ｓの処理によって出力された推定真位相マップ（φ）を示している。前述したように、この推定真位相マップ（φ）の左上には、比較的広い不連続領域が存在する。ステップＳＴ３００の処理では、コスト壁検出のために、まず、ベーシックＴＲＷ－Ｓの処理によって出力された推定真位相マップ（φ）からコストマップを生成する。具体的には、推定真位相マップ（φ）に対して、（式２１）に示す演算処理を、特定の１つの方向及び特定の１つの向きに走査しながら行うことによりコストマップを生成する。

図９（ｂ）に示す例では、推定真位相マップのＸ方向に互いに隣接する２つの画素を、Ｘ方向の上から下の向きに走査させながら（式２１）の演算処理を行うことにより、Ｘ方向（順方向）のコストマップを生成している。

その後、ステップＳＴ３００では、生成したコストマップにおいて、所定の閾値を超えるコスト画素の値が所定の範囲に連続するとき、コストマップ内にコスト壁が検出されたと判定し、該当する画素位置にコスト壁を生成する。推定真位相マップの不連続領域の外周では、位相の値が大きく変化するため、（式２１）から理解できるように、コスト画素は大きな値を示す。したがって、コスト壁が検出される位置は、主に、不連続領域の外周となる。なお、コスト画素の値が所定の閾値を超えた場合でも、その位置がコストマップ内で孤立しているような場合には、ノイズの可能性もあるため、コスト壁検出の対象からは除去している。

次のステップＳＴ３０１とステップＳＴ３０２では、コスト壁を基点とし、コストマップを生成したときの走査方向の反対側に向けてリージョングローイング処理を行うことにより、次に行われる置換処理のための誤推定領域を検出する。

図１０は、本実施形態のリージョングローイング処理の概念を説明する図である。リージョングローイング処理に先立って、先ずステップＳＴ３０１において、推定真位相マップ（φ）内の位相φ_Ａと位相φ_Ｂを入れ替えて、図１０（ａ）に例示するような反転マップを生成する。

次に、反転マップ内に、リージョングローイング処理のための基点（即ち、スタート位置）を設定する。例えば、コスト壁の中の最も値の大きな画素の位置が基点として設定される。

その後、リージョングローイング処理により、基点の画素に類似する領域を拡大していく。本実施形態のリージョングローイング処理は、領域拡大の方向に拘束条件を課した状態で行われる。例えば、コストマップが、Ｘ方向の順方向に走査することで生成された場合には、コスト壁内の位置を基点として、走査方向と反対の方向（即ち、Ｘ方向の逆方向）に向かって拡大させるように、リージョングローイング処理を行う。

また、本実施形態のリージョングローイング処理では、例えば、拡大された領域と、その領域に接する外側の領域との位相を比較し、位相の差が最も小さな外側領域の画素を拡大領域内に順次取り込むようにして、拡大領域の面積を広げていく。そして、拡大された領域と、その領域に接する外側の領域との位相差が、互いに接するすべての画素位置において所定の値以上となった場合、或いは、拡大された領域の面積が所定の面積を超えた場合に、リージョングローイング処理を停止するものとしている。

リージョングローイング処理によって拡大した領域（以下、グローイング処理後領域と呼ぶ）は、ベーシックＴＲＷ－Ｓ処理によって誤って推定した位相の領域、即ち、誤推定領域であると考えることができる。

そこで、ステップＳＴ３０３では、ベーシックＴＲＷ－Ｓ処理で用いた位相マップ（φ_Ａ）又は位相マップ（φ_Ｂ）のどちらか一方を、誤推定領域の位相で置換する。

図１１は、ステップＳＴ３０３の処理の概念を説明する図である。図１１の右図は、ステップＳＴ３０２で求めた誤推定領域を示す図であり、図１０（ｂ）と同じ図である。誤推定領域は、ベーシックＴＲＷ－Ｓ処理によって誤ったと推定される領域を示しているが、誤推定領域の内部の位相は、誤りが訂正された後の位相を示している。

したがって、図１１に例示するように、誤推定領域の内部が位相φ_Ａの場合は、オリジナルの位相マップ（φ_Ｂ）の方に対して、誤推定領域に該当する領域を、位相φ_Ｂから位相φ_Ａに置換する処理を行う。このようにして、ベーシックＴＲＷ－Ｓ処理に使用したオリジナルの位相マップ（φ_Ａ）と位相マップ（φ_Ｂ）を、位相マップ（φ_Ａ）と置換後の位相マップ（φ_Ｂ）に更新する。

ステップＳＴ３０４では、更新後の位相マップ（φ_Ａ）と位相マップ（φ_Ｂ）に対して、置換した領域の位相を固定した条件の下で、変形ＴＲＷ－Ｓアルゴリズムを実行する。置換した領域の位相を固定するという条件が付く点において、オリジナルＴＲＷ－Ｓアルゴリズムと異なっているため、変形ＴＲＷ－Ｓアルゴリズムと呼ぶものとしている。
そして、ステップＳＴ３０４では、変形ＴＲＷ－Ｓアルゴリズムの実行によって、新たな推定真位相マップ（φ）を生成し、更に、この推定真位相マップ（φ）に基づいて、新たなコストマップ（Ｃ）を生成する。

ステップＳＴ３０５では、新たなコストマップ（Ｃ）の各画素値を積算して、新たなコスト値Ｃを算出する。

図１２は、ステップＳＴ３０４、３０５の処理を説明する図である。更新後の位相マップ（φ_Ａ）と位相マップ（φ_Ｂ）に対して、上述した変形ＴＲＷ－Ｓアルゴリズムを実行することにより、オリジナルＴＲＷ－Ｓアルゴリズムで求めた推定真位相マップ（φ）（図９（ａ）参照）に対して、不連続領域が大幅に低減された更新後の推定真位相マップ（φ）（図１２の右上図）を生成することができる。

また、更新後の推定真位相マップ（φ）から算出したコストマップ（図１２の右下図）も、オリジナルＴＲＷ－Ｓアルゴリズムで算出されたコストマップに対してコスト画素値の高い領域が低減されている。このため、コストマップの全てのコスト画素値を積算して得られる新たなコスト値Ｃ２も、オリジナルＴＲＷ－Ｓアルゴリズムで算出されたコスト値Ｃ１よりも小さな値となる。

ステップＳＴ３０５の処理後は、図７のステップＳＴ２０４に戻り、全ての処理方向の順逆それぞれに対して、ステップＳＴ２０３の処理が行われたか否か（即ち、６回の処理が行われたか否か）が判定される。ＹＥＳの判定の場合には、ステップＳＴ２０５ヘ進み、オリジナルＴＲＷ－Ｓ処理で算出されたコスト値と、ステップＳＴ２０３の変形ＴＲＷ－Ｓ処理で算出された６つのコスト値が比較され、最小のコスト値に対応する推定真位相マップ（φ）が選択される。

このようにして最終的に選択された推定真位相マップ（φ）は、図５のステップＳＴ９０の処理に提供され、この推定真位相マップ（φ）を用いて、水画像と脂肪画像とが生成される。

また、必要に応じてステップＳＴ２０６の判定が行われる。ステップＳＴ２０６では、推定真位相マップ（φ）を更新したことにより、全ての処理方向に対してコスト壁が解消されたか否かの判定が行われる。そして、コスト壁が残っていると判定された場合には、上述したステップＳＴ２０１からす２０５までの処理を繰り返す。

（画像処理装置）
図１３は、実施形態に係る画像処理装置６００の構成例を示す図である。画像処理装置６００は、前述した磁気共鳴イメージング装置１のコンソール４００の機能の大部分を、磁気共鳴イメージング装置１とは独立した情報処理装置、例えば、ワークステーションやパーソナルコンピュータ、或いは、タブレット型コンピュータで行う装置である。

前述した第１画像、及び、第２画像は、磁気共鳴イメージング装置１で生成される。そして、画像処理装置６００では、磁気共鳴イメージング装置１で生成された第１画像及び第２画像を、図１３に示す外部インターフェース機能４０９、画像取得機能４１０により、任意の記録媒体や任意の通信回線を介して取得する。画像処理装置６００は、磁気共鳴イメージング装置１のコンソール４００と同様に処理回路４０を有している。

処理回路４０は、上記の画像取得機能４１０の他、第１位相算出機能４０５、第２位相算出機能４０６、改善型ＴＲＷ－Ｓ機能４０７、及び、画像生成機能４０８を実現する。これらの各機能はここまでに説明してきた各機能と同じものである。また、これらの各機能は、例えば、処理回路４０が具備するプロセッサが、記憶回路４１に記憶される所定のプログラムを実行することによって実現される。

このような画像処理装置６００によっても、前述した磁気共鳴イメージング装置１と同様の効果を得ることができる。

なお、各実施形態における再構成機能は、特許請求の範囲の記載における再構成部の一例である。また、各実施形態の説明における第１位相算出機能及び第２位相算出機能は、特許請求の範囲の記載における位相算出部の一例である。また、各実施形態の説明における改善型ＴＲＷ－Ｓ処理機能は、特許請求の範囲の記載におけるＴＲＷ－Ｓ処理部の一例である。また、各実施形態の説明における画像生成機能は、特許請求の範囲の記載における画像生成部の一例である。また、各実施形態の説明における外部インターフェース機能及び画像取得機能は、特許請求の範囲の記載における取得部の一例である

以上説明した少なくとも１つの実施形態によれば、ＴＲＷ－Ｓアルゴリズムを用いた水画像と脂肪画像の分離技術において、水成分の領域と脂肪成分の領域とを、誤りが少なく、より正確に推定できる。

本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。

１ 磁気共鳴イメージング装置
３４ シーケンスコントローラ
４０ 処理回路
４１ 記憶回路
４２ ディスプレイ
４３ 入力デバイス
４０１ 撮像条件設定機能
４０２ 再構成機能
４０３ 第１画像生成機能
４０４ 第２画像生成機能
４０５ 第１位相算出機能
４０６ 第２位相算出機能
４０７ 改善型ＴＲＷ－Ｓ処理機能
４０８ 画像生成機能
４０９ 外部インターフェース機能
４１０ 画像取得機能
６００ 画像処理装置

Claims (7)

1. 第１のエコー時間で収集されるデータから第１画像を再構成し、第２のエコー時間で収集されるデータから第２画像を再構成する再構成部と、
   前記第１画像と前記第２画像とから、画素ごとに、一方が真であり他方が偽である第１位相及び第２位相を算出し、第１位相マップ及び第２位相マップを夫々生成する位相算出部と、
   前記第１位相マップ及び第２位相マップから、ＴＲＷ－Ｓ（sequential tree-reweighted message-passing）アルゴリズムを用いて、被検体の撮像時の静磁場分布に対応する真の位相分布を推定した推定真位相マップを生成するＴＲＷ－Ｓ処理部と、
   前記第１及び第２画像と推定真位相マップとを用いて、水画像と脂肪画像とを生成する画像生成部と、
   を備え、
   前記ＴＲＷ－Ｓ処理部は、
   （ａ）前記推定真位相マップの誤推定領域を検出し、
   （ｂ）検出した前記誤推定領域が前記第１位相であるべき領域を前記第２位相の領域であると誤って推定している場合には、第２位相マップの該当領域を前記第１位相に置換し、検出した前記誤推定領域が前記第２位相であるべき領域を前記第１位相の領域であると誤って推定している場合には、前記第１位相マップの該当領域を前記第２位相に置換し、
   （ｃ）置換した領域の位相を固定した状態で、前記ＴＲＷ－Ｓアルゴリズムを実行して、前記推定真位相マップを更新し、
   前記画像生成部は、更新された前記推定真位相マップと、前記第１及び第２画像とを用いて、前記水画像と前記脂肪画像とを生成する、
   磁気共鳴イメージング装置。
2. 前記ＴＲＷ－Ｓ処理部は、
   前記推定真位相マップの隣接する位相の滑らかさを表す指標値を夫々の画素値とするコストマップを生成し、
   前記コストマップの画素値の値に基づいて、前記誤推定領域を検出する、
   請求項１に記載の磁気共鳴イメージング装置。
3. 前記ＴＲＷ－Ｓ処理部は、
   前記コストマップを、初期の推定真位相マップと、更新後の推定真位相マップの夫々に対して生成し、更に、前記コストマップの各画素値を積算したコスト値を、前記初期の更推定真位相マップと、前記更新後の推定真位相マップの夫々に対して算出し、
   前記更新後の推定真位相マップのコスト値が、前記初期の推定真位相マップのコスト値よりも小さくなったとき、前記更新後の推定真位相マップを最終的な推定真位相マップとする、
   請求項２に記載の磁気共鳴イメージング装置。
4. 前記ＴＲＷ－Ｓ処理部は、
   （ａ）前記推定真位相マップを所定の１つの方向に沿って走査することによって前記コストマップを生成し、
   （ｂ）前記コストマップにおいて、所定の閾値を超える画素値が所定の範囲に連続するとき、前記コストマップ内にコスト壁を生成し、
   （ｃ）前記コスト壁を基点とし、前記走査方向の反対側に向けてリージョングローイング処理を行うことにより、前記誤推定領域を検出し、
   （ｄ）前記誤推定領域の置換によって前記推定真位相マップを更新する、
   請求項２に記載の磁気共鳴イメージング装置。
5. 前記ＴＲＷ－Ｓ処理部は、
   （ａ）初期の推定真位相マップに対するコストマップを生成すると共に、前記初期の推定真位相マップを、互いに直交する３軸の夫々に対して順方向及び逆方向に走査することによって、６つの更新後の推定真位相マップと６つのコストマップを生成し、
   （ｂ）前記初期の推定真位相マップのコストマップと、前記６つのコストマップに対して、各コストマップの各画素値を積算したコスト値を夫々算出し、
   （ｃ）コスト値が最も小さいコストマップに対応する推定真位相マップを、最終的な推定真位相マップとして出力する、
   請求項４に記載の磁気共鳴イメージング装置。
6. 第１のエコー時間で収集されるデータから再構成された第１画像と、第２のエコー時間で収集されるデータから再構成された第２画像とを取得する画像取得部と、
   前記第１画像と前記第２画像とから、画素ごとに、一方が真であり他方が偽である第１位相及び第２位相を算出し、第１位相マップ及び第２位相マップを夫々生成する位相算出部と、
   前記第１位相マップ及び第２位相マップから、ＴＲＷ－Ｓ（sequential tree-reweighted message-passing）アルゴリズムを用いて、被検体の撮像時の静磁場分布に対応する真の位相分布を推定した推定真位相マップを生成するＴＲＷ－Ｓ処理部と、
   前記第１及び第２画像と推定真位相マップとを用いて、水画像と脂肪画像とを生成する画像生成部と、
   を備え、
   前記ＴＲＷ－Ｓ処理部は、
   （ａ）前記推定真位相マップの誤推定領域を検出し、
   （ｂ）検出した前記誤推定領域が前記第１位相であるべき領域を前記第２位相の領域であると誤って推定している場合には、第２位相マップの該当領域を前記第１位相に置換し、検出した前記誤推定領域が前記第２位相であるべき領域を前記第１位相の領域であると誤って推定している場合には、前記第１位相マップの該当領域を前記第２位相に置換し、
   （ｃ）置換した領域の位相を固定した状態で、前記ＴＲＷ－Ｓアルゴリズムを実行して、前記推定真位相マップを更新し、
   前記画像生成部は、更新された前記推定真位相マップと、前記第１及び第２画像とを用いて、前記水画像と前記脂肪画像とを生成する、
   画像処理装置。
7. 第１のエコー時間で収集されるデータから再構成された第１画像と、第２のエコー時間で収集されるデータから再構成された第２画像とを取得し、
   前記第１画像と前記第２画像とから、画素ごとに、一方が真であり他方が偽である第１位相及び第２位相を算出し、第１位相マップ及び第２位相マップを夫々生成し、
   前記第１位相マップ及び第２位相マップから、ＴＲＷ－Ｓ（sequential tree-reweighted message-passing）アルゴリズムを用いて、被検体の撮像時の静磁場分布に対応する真の位相分布を推定した推定真位相マップを生成し、
   前記第１及び第２画像と推定真位相マップとを用いて、水画像と脂肪画像とを生成する、ステップを備え、
   前記推定真位相マップを生成するステップでは、
   （ａ）前記推定真位相マップの誤推定領域を検出し、
   （ｂ）検出した前記誤推定領域が前記第１位相であるべき領域を前記第２位相の領域であると誤って推定している場合には、第２位相マップの該当領域を前記第１位相に置換し、検出した前記誤推定領域が前記第２位相であるべき領域を前記第１位相の領域であると誤って推定している場合には、前記第１位相マップの該当領域を前記第２位相に置換し、
   （ｃ）置換した領域の位相を固定した状態で、前記ＴＲＷ－Ｓアルゴリズムを実行して、前記推定真位相マップを更新し、
   前記水画像と脂肪画像とを生成するステップでは、更新された前記推定真位相マップと、前記第１及び第２画像とを用いて、前記水画像と前記脂肪画像とを生成する、
   画像処理方法。

Images