JP2022002762A - 磁気共鳴イメージングで得た画像の処理方法、画像処理プロブラム、及び、計算機 - Google Patents

Abstract

【課題】脳表領域を含めた脳全体の局所位相分布を高精度に算出する磁気共鳴イメージング装置を提供する。【解決手段】磁気共鳴イメージング装置は、被検体が配置される空間に静磁場を生成する静磁場生成磁石と、前記被検体に高周波磁場を送信する送信部と、前記高周波磁場の送信により前記被検体から発生する核磁気共鳴信号を受信する受信部と、前記核磁気共鳴信号に位置情報を付加する傾斜磁場を印加する傾斜磁場印加部と、前記送信部、前記受信部、及び前記傾斜磁場印加部を制御するとともに、前記核磁気共鳴信号に対して演算処理を行う計算機と、前記計算機により演算処理された画像を表示する表示部と、を備え、前記計算機は、所定の演算処理を行う。【選択図】図３

Description

本発明は、磁気共鳴イメージング（MRI:magnetic resonance imaging）技術に関する。特に、Gradient Echo法で撮像された複素画像を用いて、生体組織間の磁化率差に起因した局所位相分布を算出する画像処理技術に関する。

磁気共鳴イメージング装置は、静磁場内に配置される被検体の水素原子核（プロトン）が、特定周波数の高周波磁場に共鳴する核磁気共鳴（NMR:nuclear magnetic resonance）現象を利用した非侵襲な医用画像診断装置である。被検体から発生する核磁気共鳴信号は、プロトン密度や緩和時間など様々な物性値によって変化するため、磁気共鳴イメージング装置は、核磁気共鳴信号に基づいて、生体組織の構造や組成、細胞性状などの様々な生体情報を描出できる。

近年、磁気共鳴イメージング装置により測定可能な物性値として、生体組織間の磁化率差が注目されている。磁化率とは、静磁場中の物質が磁気分極（磁化）する度合いを表した物性値である。生体内には、静脈血中のデオキシヘモグロビンや鉄たんぱく質などの常磁性体と、生体組織の大部分を占める水や石灰化の基となるカルシウムなどの反磁性体とが存在する。磁気共鳴イメージング装置は、生体組織間の磁化率差を画像化することができるため、脳虚血疾患の診断やがんの放射線治療効果の予測、神経変性疾患の鑑別、などへの可能性が期待されている。生体組織間の磁化率差を画像化する方法としては、例えば、定量的磁化率マッピング（QSM:quantitative susceptibility mapping）法、磁化率強調画像化（SWI:Susceptibility Weighted Imaging）法、などが知られている。

磁気共鳴イメージング装置は、生体組織間の磁化率差に起因する局所位相分布（局所的な位相変化の分布）を算出し、当該局所位相分布を用いてＱＳＭ法により定量的磁化率分布を算出する、或いは、当該局所位相分布を用いてＳＷＩ法により磁化率強調画像を算出する。磁気共鳴イメージング装置は、当該局所位相分布を算出する際、位相折り返し補正処理を行って総位相分布を算出し、総位相分布から被検体の形状などに起因して生じる背景位相分布を除去する背景位相除去処理を行う。背景位相除去処理における代表的な方法として、以下の二つの方法が提案されている。

一つ目の方法は、ＰＤＦ（Projection onto Dipole Fieled）法と呼ばれる方法である（例えば、特許文献１参照）。ＰＤＦ法は、計測した脳内部の位相分布から脳外部の磁化率分布を推定し、該磁化率分布によって生じる脳内部の位相分布を背景位相分布として算出し、総位相分布から背景位相分布を除去することで局所位相分布を算出する。
二つ目の方法は、ＳＨＡＲＰ（Sophisticated Harmonic Artifact Reduction for Phase date）法と呼ばれる方法である（例えば、非特許文献１参照）。ＳＨＡＲＰ法は、背景位相分布が球面調和関数近似できることを利用する。即ち、任意の点における値が、その点を中心とした任意の半径を有する球（カーネル）内の点における値の平均に等しい、という平均値の性質（Spherical mean value property）を用いて、総位相分布から背景位相分布と局所位相分布とを分離する。

米国特許第９４４８２８９号明細書

Ｆｅｒｄｉｎａｎｄ Ｓｃｈｗｅｓｅｒ他、 "Ｑｕａｎｔｉｔａｔｉｖｅ ｉｍａｇｉｎｇ ｏｆ ｉｎｔｒｉｎｓｉｃ ｍａｇｎｅｔｉｃ ｔｉｓｓｕｅ ｐｒｏｐｅｒｔｉｅｓ ｕｓｉｎｇ ＭＲＩ ｓｉｇｎａｌ ｐｈａｓｅ： Ａｎ ａｐｐｒｏａｃｈ ｔｏ ｉｎ ｖｉｖｏ ｂｒａｉｎ ｉｒｏｎ ｍｅｔａｂｏｌｉｓｍ？"、 Ｎｅｕｒｏｉｍａｇｅ、２０１１年、５４巻、２７８９−２８０７頁

しかしながら、特許文献１に記載の方法では、脳外部の推定磁化率で再現できない部分（例えば、副鼻腔など）の背景位相分布を除去することができないという問題がある。また、非特許文献１に記載の方法では、脳内部と脳外部との境界部分（脳表領域）における局所位相分布の算出精度が低下してしまうという問題がある。

本発明は、上記した課題を解決するためになされたものであり、脳表領域を含めた脳全体の局所位相分布を高精度に算出する磁気共鳴イメージング装置を提供することを課題とする。

前記課題を解決するために、本開示の実施形態に係る磁気共鳴イメージング装置は、被検体が配置される空間に静磁場を生成する静磁場生成磁石と、前記被検体に高周波磁場を送信する送信部と、前記高周波磁場の送信により前記被検体から発生する核磁気共鳴信号を受信する受信部と、前記核磁気共鳴信号に位置情報を付加する傾斜磁場を印加する傾斜磁場印加部と、前記送信部、前記受信部、及び前記傾斜磁場印加部を制御するとともに、前記核磁気共鳴信号に対して演算処理を行う計算機と、前記計算機により演算処理された画像を表示する表示部と、を備え、前記計算機は、前記核磁気共鳴信号によって構成された複素画像から前記被検体の対象組織を抽出した対象組織抽出画像を算出する対象組織抽出部と、前記複素画像から位相分布を算出する位相分布算出部と、前記位相分布の位相折り返しを補正する位相折り返し補正部と、前記対象組織抽出画像に基づいて、補正後の位相分布を、前記対象組織の所定領域と前記所定領域以外の領域とに分割する領域分割部と、前記所定領域の背景位相分布を算出する背景位相分布算出部と、前記補正後の位相分布と前記背景位相分布とを組み合わせて、前記対象組織の所定領域における局所位相分布を算出する局所位相分布算出部と、を備えることを特徴とする。

脳表領域を含めた脳全体の局所位相分布を高精度に算出する磁気共鳴イメージング装置を提供できる。

（ａ）は、本実施形態に係る垂直磁場方式の磁気共鳴イメージング装置の外観図、（ｂ）は、本実施形態に係る水平磁場方式の磁気共鳴イメージング装置の外観図、（ｃ）は、本実施形態に係る開放感を高めた磁気共鳴イメージング装置の外観図である。 本実施形態に係る磁気共鳴イメージング装置の概略構成を示すブロック図である。 本実施形態に係る計算機の機能ブロック図である。 本実施形態に係る計算機における演算処理の流れを説明するフローチャートである。 ＲＳＳＧ（RF-spoiled-Steady-state Acquisition with Rewound Gradient-Echo）シーケンスのパルスシーケンスを説明するための図である。 本実施形態に係る局所位相分布算出部における局所位相分布算出処理の流れを説明するフローチャートである。 （ａ）は、本実施形態に係る全領域局所位相分布算出部で用いる重み画像Ｗ_１を説明するための図、（ｂ）は、本実施形態に係る全領域局所位相分布算出部で用いる重み画像Ｗ_２を説明するための図である。 （ａ）は、本実施形態に係る計算機によって演算処理が施された画像を示す図、（ｂ）は、（ａ）の一部領域における拡大図、（ｃ）は、従来の計算機によって演算処理が施された画像を示す図、（ｄ）は、（ｃ）の一部領域における拡大図である。

本発明を適用する実施形態の一例を説明する。以下、実施形態を説明するための全図において、特に断らない限り、同一機能を有するものは同一符号を付し、繰り返しの説明は省略する。また、以下の記述により本発明が限定されるものではない。

≪ＭＲＩ装置の全体構成≫
図１は、本実施形態に係るＭＲＩ装置の外観図である。図１（ａ）は、ソレノイドコイルで静磁場を生成するトンネル型磁石を用いた水平磁場方式のＭＲＩ装置１００である。図１（ｂ）は、開放感を高めるために磁石を上下に分離したオープン型の垂直磁場方式のＭＲＩ装置１２０である。また、図１（ｃ）は、図１（ａ）と同じトンネル型磁石を用い、磁石の奥行を短くし、かつ、斜めに傾けることによって、開放感を高めたＭＲＩ装置１３０である。なお、本実施形態において、ＭＲＩ装置の構成は、特に限定されるものではなく、図１に示す何れの構成であってもよいし、公知の構成であってもよい。

以下、特に区別する必要がない場合は、ＭＲＩ装置１００で代表する。ＭＲＩ装置１００において、ＭＲＩ装置１００の静磁場方向をｚ方向、ｚ方向に垂直な２方向のうち、測定対称である被検体１０１を載置するベッド面に平行な方向をｘ方向、他方向をｙ方向とする。

図２は、本実施形態に係るＭＲＩ装置の概略構成を示すブロック図である。図２に示すように、ＭＲＩ装置１００は、被検体１０１が配置される空間に静磁場を生成する静磁場生成磁石１０２と、被検体１０１に高周波磁場を送信する送信用高周波コイル（送信部）１０５と、高周波磁場の送信により被検体１０１から発生する核磁気共鳴信号を受信する受信用高周波コイル（受信部）１０６と、核磁気共鳴信号に位置情報を付加する傾斜磁場を印加する傾斜磁場印加部１０３と、送信部１０５、受信部１０６、及び傾斜磁場印加部１０３を制御するとともに、核磁気共鳴信号に対して演算処理を行う計算機１０９と、計算機１０９により演算処理された画像を表示する表示部１１０と、を備える。
更に、ＭＲＩ装置１００は、シムコイル１０４、送信機１０７、受信機１０８、外部記憶装置１１１、傾斜磁場用電源部１１２、シム用電源部１１３、シーケンス制御装置１１４、入力装置１１５、などを備える。

静磁場生成磁石１０２は、被検体１０１が配置される空間に静磁場を生成する磁石であり、ＭＲＩ装置の構成に応じて、種々の形態を採用することができる。傾斜磁場印加部１０３は、核磁気共鳴信号に空間的な位置情報を付加するために、被検体１０１が配置される空間において、ｘ方向、ｙ方向、ｚ方向の各方向に傾斜磁場を印加するコイルである。各方向に印加される傾斜磁場（ｘ方向傾斜磁場、ｙ方向傾斜磁場、ｚ方向傾斜磁場）は、例えば、位相エンコード傾斜磁場、リードアウト傾斜磁場、スライス傾斜磁場である。

シムコイル１０４は、被検体１０１が配置される空間に生成される静磁場の不均一性を打ち消して除去するなど、静磁場分布を調整するために使用される電流コイルである。送信部１０５は、被検体１０１の水素原子核を励起させる高周波磁場を、被検体１０１に送信する送信用高周波コイルである。受信部１０６は、高周波磁場の送信により被検体１０１から発生する核磁気共鳴信号を受信する受信用高周波コイルである。なお、送信部１０５と受信部１０６とは、別々の構成としてもよいし、兼用する構成としてもよい。

送信機１０７は、シーケンス制御装置１１４によって制御され、送信部１０５へと送信する高周波磁場を生成する。受信機１０８は、シーケンス制御装置１１４によって制御され、受信部１０６によって受信された核磁気共鳴信号を検波し、複素信号を計算機１０９へと送信する。なお、送信部１０５と受信部１０６とは、別々の構成としてもよいし、兼用する構成としてもよい。

計算機１０９は、ＣＰＵ（Central Processing Unit）、メモリ、記憶装置などを備える情報処理装置であり、計算機１０９には、表示部１１０、外部記憶装置１１１、入力装置１１５などが接続される。計算機１０９は、記憶装置に記憶される制御プログラムを読み出して、ワークエリアに展開し、当該制御プログラムを実行することで、ＭＲＩ装置１００全体の動作を制御するとともに、各種演算処理を行う。なお、計算機１０９は、ＭＲＩ装置１００の一部であってもよいし、ＭＲＩ装置１００から独立した情報処理装置であって、ＭＲＩ装置１００とデータの送受信が可能な情報処理装置であってもよい。更に、計算機１０９が実現する各種機能のうち、全部または一部の機能は、ＡＳＩＣ（Applicat
ion Specific Integrated Circuit）、ＦＰＧＡ（field-programmable gate array）などのハードウェアによって実現されていてもよい。

表示部１１０は、例えば、液晶ディスプレイ、有機ＥＬディスプレイ、等で構成される。表示部１１０は、計算機１０９から取得したデータに基づいて、ディスプレイ画面上に、演算結果、各種画像などを表示する。表示部１１０は、例えば、後述する局所位相分布、定量的磁化率分布、磁化率強調画像などをディスプレイ画面上に表示する。

外部記憶装置１１１は、記憶装置とともに、計算機１０９が行う各種演算処理に使用されるデータ、演算処理中に生成されるデータ、各種演算処理により取得されるデータ、入力装置１１５を介して入力される各種条件、各種パラメータ、などを記憶する。

傾斜磁場用電源部１１２は、傾斜磁場印加部１０３の駆動用電源であり、シーケンス制御装置１１４からの制御指示に基づいて、傾斜磁場印加部１０３の駆動を制御する。シム用電源部１１３は、シムコイル１０４の駆動用電源であり、シーケンス制御装置１１４からの制御指示に基づいて、シムコイル１０４の駆動を制御する。

シーケンス制御装置１１４は、送信機１０７、受信機１０８、傾斜磁場用電源部１１２、シム用電源部１１３、などＭＲＩ装置１００の各部の動作を制御する。例えば、シーケンス制御装置１１４は、外部記憶装置１１１や記憶装置に記憶される設定情報（例えば、後述のパルスシーケンスと呼ばれるタイムチャート）に基づいて、高周波磁場の印加タイミング、傾斜磁場（ｘ方向傾斜磁場、ｙ方向傾斜磁場、ｚ方向傾斜磁場）の印加タイミング、核磁気共鳴信号の受信タイミング、などを制御する。

入力装置１１５は、操作者による入力操作を受け付ける機能を備えており、例えば、キーボード、カメラ、マウス、タッチパネル、音声入力を受け付けるマイク、などで構成される。入力装置１１５は、計算機１０９が行う各種演算処理や計測処理などに必要な各種条件、各種パラメータ、などを操作者が入力するためのインタフェースである。なお、各種条件、各種パラメータ、などは、予め設定しておくことも可能であるし、操作者が入力装置１１５を介して適宜設定することも可能である。

≪計算機の構成≫
次に、図３を参照して、本実施形態に係る計算機１０９の構成について説明する。本実施形態では対象組織として被検体１０１の脳を適用する場合を一例に挙げて説明する。

ここで、計算機１０９の構成の説明に先立ち、ＳＨＡＲＰ法を用いた一般的な背景位相除去処理の概要について説明する。
まず、ＭＲＩ装置は、Gradient Echo法を用いて、任意のエコー時間（高周波磁場を送信してから核磁気共鳴信号を受信するまでの時間）で撮像された複素画像から位相分布Ｐを算出する。そして、ＭＲＩ装置は、位相分布Ｐに生じている−π〜+πの範囲に折り返された位相を補正する位相折り返し補正処理を行う。
位相折り返し補正処理後の位相分布（総位相分布）Ｐ_{ｔｏｔａｌ}は、対象組織の形状などに起因する背景位相分布Ｐ_ｂｋｇｒ、生体組織間の磁化率差に起因する局所位相分布Ｐ_{ｌｏｃａｌ}を用いて、式（１）で表される。

そして、ＭＲＩ装置は、任意画素における画素値と、その画素を中心とした所定の半径を有する球(カーネル)内の画素における画素値の平均とが略等しいという平均値の性質を用いて、総位相分布Ｐ_{ｔｏｔａｌ}から、背景位相分布Ｐ_ｂｋｇｒと局所位相分布Ｐ_{ｌｏｃａｌ}とを分離する。
背景位相分布Ｐ_ｂｋｇｒは、平均値の性質によって、式（２）を満たす。

δはデルタ関数、ρ_ｒは半径ｒの球状カーネル、Ａは畳み込み積分をそれぞれ表す。また、Ｍ_{ｓｈｒｉｎｋ}は、被検体１０１の対象組織を抽出した対象組織抽出画像Ｍから、球状カーネルρ_ｒの半径ｒに応じて縮小されたマスクを表す。
式（２）で表される畳み込み積分Ａは、フーリエ変換行列Ｆ、逆フーリエ変換行列Ｆ^−１を用いることで、式（３）に変換できる。

但し、Ｃ＝Ｆ(δ−ρ_ｒ)
式（３）に式（１）を代入することで、式（４）の関係式が得られる。

そして、ＭＲＩ装置は、総位相分布Ｐ_{ｔｏｔａｌ}から、式（４）の関係式を満たす局所位相分布Ｐ_{ｌｏｃａｌ}を算出し、最小二乗処理によって、式（５）で表される局所位相分布Ｐ’_{ｌｏｃａｌ}を算出する。

即ち、ＳＨＡＲＰ法を用いた一般的な背景位相除去処理では、マスクＭ_{ｓｈｒｉｎｋ}の範囲における局所位相分布しか算出することができず、脳内部と脳外部との境界部分（脳表領域）における局所位相分布の算出精度が低下するという問題が生じていた。

そこで、本実施形態に係るＭＲＩ装置１００は、総位相分布を、脳表領域と脳内部の領域とに領域分割し、総位相分布と脳表領域における背景位相分布とを組み合わせることで、脳表領域を含めた脳全体の局所位相分布を算出する。これにより、脳表領域における局所位相分布の算出精度を向上させることができる。以下、これを実現する計算機１０９の機能構成について説明する。

計算機１０９は、計測制御部３１０と、画像再構成部３２０と、局所位相分布算出部３３０と、所望分布算出部３４０と、を備える。更に、局所位相分布算出部３３０は、対象組織抽出部３３１、位相分布算出部３３２、位相折り返し補正部３３３、領域分割部３３４、背景位相分布算出部３３５、全領域局所位相分布算出部３３６と、を備える。

計測制御部３１０は、例えば、Gradient Echo法を用いて、送信部１０５から被検体１０１への高周波磁場の送信により、被検体１０１から発生する任意のエコー時間の核磁気共鳴信号を複素信号として計測する。計測制御部３１０は、予め設定されたパルスシーケンスに従って、シーケンス制御装置１１４を制御することで、高周波磁場の印加タイミング、傾斜磁場の印加タイミング、核磁気共鳴信号の受信タイミング、などを制御する。なお、計測制御部３１０は、血流などの流れの影響を補償する流速補正（Flow Compensation）傾斜磁場パルスを、ｘ方向、ｙ方向、ｚ方向の各方向に印加してもよい。

画像再構成部３２０は、計測制御部３１０によって計測された複素信号（核磁気共鳴信号）を、例えば、ｋｒ軸、ｋｐ軸、ｋｓ軸を座標軸とする３次元のｋ空間に配置してフーリエ変換することにより、任意画素の画素値が複素数で表される複素画像を再構成する。画像再構成部３２０は、ｋ空間に核磁気共鳴信号を充填する際、ｋ空間の座標軸（例えば、ｋｒ軸）に平行に、１行ずつ核磁気共鳴信号を充填していく方法を用いてもよいし（Cartesian充填法）、ｋ空間の原点を中心に、角度を変化させながら放射状に核磁気共鳴信号を充填していく方法を用いてもよい（Radial Scan充填法）。

局所位相分布算出部３３０は、画像再構成部３２０によって再構成された複素画像に基づいて、対象組織の全領域における局所位相分布を算出する。具体的には、対象組織抽出部３３１が、複素画像から被検体１０１の対象組織を抽出した対象組織抽出画像Ｍを生成する。そして、位相分布算出部３３２が、複素画像から位相分布Ｐを算出する。そして、位相折り返し補正部３３３が、位相分布Ｐの位相折り返しを補正し、総位相分布Ｐ_{ｔｏｔａｌ}を算出する。そして、領域分割部３３４が、対象組織抽出画像Ｍに基づいて、総位相分布Ｐ_{ｔｏｔａｌ}を少なくとも２つ以上の領域（例えば、所定領域、所定領域以外の領域）に分割する。そして、背景位相分布算出部３３５が、所定領域における背景位相分布Ｐｅｄｇｅを算出する。そして、全領域局所位相分布算出部３３６が、総位相分布Ｐ_{ｔｏｔａｌ}と、所定領域における背景位相分布Ｐ_ｅｄｇｅとを組み合わせることで、所定領域を含めた対象組織の全領域における局所位相分布Ｐ’_{ｌｏｃａｌ}を算出する。

所望分布算出部３４０は、局所位相分布算出部３３０によって算出された局所位相分布Ｐ’_{ｌｏｃａｌ}を用いて、定量的磁化率分布や磁化率強調画像などの所望分布を算出し、当該所望分布を表示部１１０へと出力する。

上述のように、本実施形態に係る計算機１０９は、複素画像から対象組織のみを抽出し、総位相分布を対象組織の所定領域と所定領域以外の領域とに分割して、所定領域における背景位相分布を算出し、当該背景位相分布を利用して、対象組織の全領域における局所位相分布を算出する。これにより、カーネルの半径に応じて脳表領域で誤差が大きくなり、脳表領域における局所位相分布の算出精度が低下するという従来のＳＨＡＲＰ法に生じていた問題点を回避し、脳表領域を含めた脳全体における局所位相分布を高精度に算出することができる。

≪計算機における演算処理の流れ≫
次に、図４乃至図８を参照して、本実施形態に係る計算機１０９が行う演算処理の流れについて説明する。図４は、本実施形態に係る計算機１０９における演算処理の流れを説明するフローチャートである。図５は、ＲＳＳＧシーケンスのパルスシーケンスを説明するための図である。図６は、本実施形態に係る局所位相分布算出部３３０における局所位相分布算出処理の流れを説明するフローチャートである。図７は、本実施形態に係る全領域局所位相分布算出部３３６で用いる重み画像を説明するための図である。図８は、本実施形態に係る表示部１１０に表示される画像の一例を示す図である。

［計測：Ｓ１００１］
ステップＳ１００１において、計算機１０９は計測処理を行う。

操作者が、入力装置１１５を介して各種パラメータ（例えば、核磁気共鳴信号の数、エコー時間、エコー時間の間隔など）を設定し、計算機１０９に撮像開始指示が入力されると、計測制御部３１０は、予め設定されたパルスシーケンスに従って、シーケンス制御装置１１４を制御する。そして、計測制御部３１０は、任意のエコー時間の核磁気共鳴信号を複素信号として計測する。

ここで、図５を参照して、計測制御部３１０が計測に使用するGradient Echo型のパルスシーケンスについて説明する。図５に示すＲＳＳＧ（RF-spoiled-Steady-state Acquisition with Rewound Gradient-Echo）シーケンス５１０は、生体組織間の磁化率差によって生じる局所的な磁場変化に鋭敏なパルスシーケンスである。なお、計測制御部３１０が計測に使用するシーケンスは、ＲＳＳＧシーケンス５１０に限定されるものではない。

図５において、ＲＦは高周波磁場を示し、Ｇｓはスライス傾斜磁場を示し、Ｇｐは位相エンコード傾斜磁場を示し、Ｇｒはリードアウト傾斜磁場を示している。また、ＴＥはエコー時間を示し、ＴＲは繰り返し時間を示している。更に、５０１はスライス傾斜磁場パルス、５０２は高周波磁場パルス、５０３、５０８はスライスエンコード傾斜磁場パルス、５０４、５０９は位相エンコード傾斜磁場パルス、５０５、５０６はリードアウト傾斜磁場パルス、５０７は核磁気共鳴信号を示している。なお、ハイフン以下の数字は、計測制御部３１０が行う繰り返しの計測が何回目であるかを示しており、例えば、「−２」は、計測制御部３１０が行う繰り返しの計測が２回目であることを示している。

計測制御部３１０は、ＲＳＳＧシーケンス５１０に従って、シーケンス制御装置１１４を制御することで、スライス傾斜磁場パルス５０１、高周波磁場パルス５０２、スライスエンコード傾斜磁場パルス５０３、位相エンコード傾斜磁場パルス５０４、リードアウト傾斜磁場パルス５０５、リードアウト傾斜磁場パルス５０６、スライスエンコード傾斜磁場パルス５０８、位相エンコード傾斜磁場パルス５０９、などの印加タイミングを適切に制御することで、画像再構成に必要なデータを計測する。

図５に示すＲＳＳＧシーケンス５１０では、繰り返し時間ＴＲ内に、１つの核磁気共鳴信号５０７を計測する。
まず、静磁場内に被検体１０１が配置されると、静磁場強度に応じて被検体１０１内の水素原子核が、所定の周波数を持つため、送信部１０５は、被検体１０１の計測領域に、当該周波数と一致するような高周波磁場パルス５０２−１を送信し、被検体１０１の水素原子核を励起させる。この際、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｚ方向に、スライス傾斜磁場パルス５０１−１を印加し、被検体１０１における所定のスライス（３次元画像の撮像位置）での磁化を選択励起させる。

なお、計測制御部３１０は、繰り返しの計測を行う際、高周波磁場パルス５０２の位相を、所定角度（例えば、１１７度、１２２度）ずつ変更する。例えば、繰り返しの計測が１回目の場合、所定角度を１１７度とし、繰り返しの計測が２回目の場合、角度を１１７×２度とし、繰り返しの計測が３回目の場合、角度を１１７×３度とする。これにより、次の高周波磁場パルスを印加する際、余分な信号成分を残すことなく適切な計測を行うことが可能になる。

続いて、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｚ方向に、スライスエンコード傾斜磁場パルス５０３−１を印加して、磁化の位相に、スライスエンコード方向の位置情報を付加する。また、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｘ方向に、位相エンコード傾斜磁場パルス５０４−１を印加して、磁化の位相に、位相エンコード方向の位置情報を付加する。即ち、傾斜磁場印加部１０３は、核磁気共鳴信号の読み出しに先だって、ｚ方向及びｘ方向における磁化の位相に、位置情報を付加する。

スライスエンコード傾斜磁場パルス５０３−１は、繰り返し時間ＴＲ毎に、磁場強度が変化する。例えば、図５に示す７本の横線は、異なる強さの傾斜磁場を示しており、上方向の矢印に従って、正の傾斜磁場が強くなることを示している。同様に、位相エンコード傾斜磁場パルス５０４−１も、繰り返し時間ＴＲ毎に、磁場強度が変化する。例えば、図５に示す７本の横線は、異なる強さの傾斜磁場を示しており、上方向の矢印に従って、正の傾斜磁場が強くなることを示している。

なお、図５では、スライスエンコード傾斜磁場パルス５０３−１及び位相エンコード傾斜磁場パルス５０４−１における磁場強度が、７通りに変化する場合を一例に挙げて説明しているが、これに限定されるものではない。

続いて、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｙ方向に、ディフェーズ用のリードアウト傾斜磁場パルス５０５−１を印加して、画素内における磁化の位相を分散させる。そして、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｙ方向に、リードアウト傾斜磁場パルス５０６−１を印加して、磁化の位相に、リードアウト方向の位置情報を付加する。これにより、計測制御部３１０は、繰り返し時間ＴＲ内に、１つの核磁気共鳴信号５０７を計測することが可能になる。

最後に、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｚ方向に、リフェーズ用のスライスエンコード傾斜磁場パルス５０８−１を印加し、更に、被検体１０１の計測領域におけるｘ方向に、リフェーズ用の位相エンコード傾斜磁場パルス５０９−１を印加して、ディフェーズされた画素内における磁化の位相を収束させる。

スライスエンコード傾斜磁場パルス５０８−１は、繰り返し時間ＴＲ毎に、磁場強度が変化する。例えば、図５に示す７本の横線は、異なる強さの傾斜磁場を示しており、下方向の矢印に従って、負の傾斜磁場が強くなることを示している。同様に、位相エンコード傾斜磁場パルス５０９−１も、繰り返し時間ＴＲ毎に、磁場強度が変化する。例えば、図５に示す７本の横線は、異なる強さの傾斜磁場を示しており、下方向の矢印に従って、負の傾斜磁場が強くなることを示している。

計測制御部３１０は、上述の手順を、高周波磁場パルス５０２の位相、スライスエンコード傾斜磁場パルス５０３、５０８の強度、位相エンコード傾斜磁場パルス５０４、５０
９の強度、などを変化させながら繰り返し行い、画像再構成に必要とされる核磁気共鳴信号５０７を計測する。なお、ＲＳＳＧシーケンス５１０により得られる複素画像の絶対値成分は、エコー時間ＴＥが短く設定されるとＴ１(縦緩和時間)強調画像となり、エコー時間ＴＥが長く設定されると画素内の位相分散を反映したＴ２^＊（見かけの横緩和時間）強調画像となる。本実施形態では、エコー時間ＴＥを長く設定し、Ｔ２^＊強調画像の撮像を行う。

［画像再構成：Ｓ１００２］
ステップＳ１００２において、計算機１０９は画像再構成処理を行う。

画像再構成部３２０は、ステップＳ１００１で計測された核磁気共鳴信号を、ｋ空間上に配置してフーリエ変換することにより、複素画像を再構成する。撮像方法としては、特に限定されるものではなく、カーテシアン撮像、ノンカーテシアン撮像、マルチエコーエコープラナーイメージング法、などの公知の方法を適用することができる。

［局所位相分布算出処理：Ｓ１００３］
ステップＳ１００３において、計算機１０９は局所位相分布算出処理を行う。

ここで、図６を参照して、局所位相分布算出処理の流れについて詳細に説明する。局所位相分布算出部３３０は、ステップＳ１００２で再構成された複素画像に基づいて、対象組織の全領域における局所位相分布を算出する。

ステップＳ１１０１において、対象組織抽出部３３１は、画像再構成部３２０によって構成された複素画像の絶対値成分（絶対値画像）における画素値の大きさに基づいて、被検体１０１の対象組織を抽出した対象組織抽出画像Ｍを算出する。対象組織抽出部３３１は、例えば、判別分析法などの公知の方法を用いて、絶対値画像の画素値のヒストグラムから、ノイズ成分と信号成分とを分離し、被検体１０１の対象組織を抽出した対象組織抽出画像Ｍ（例えば、脳内部の画素値を１、脳外部の画素値を０とした画像）を算出する。なお、対象組織抽出部３３１は、例えば、領域拡大（Region Growing）法などの公知の方法を用いて、複素画像から磁化率分布の算出に不要な領域（例えば、皮下脂肪などの領域）を除去して、磁化率分布の算出に必要な領域である対象組織の実質領域のみを抽出した対象組織抽出画像Ｍを算出してもよい。

ステップＳ１１０２において、位相分布算出部３３２は、画像再構成部３２０によって再構成された複素画像の偏角成分に基づいて、位相分布Ｐを算出する。

ステップＳ１１０３において、位相折り返し補正部３３３は、例えば、領域拡大法など公知の方法を用いて、ステップＳ１１０２で算出された位相分布Ｐに生じている−π〜＋πの範囲に折り返された位相を補正し、総位相分布Ｐ_{ｔｏｔａｌ}を算出する。

ステップＳ１１０４において、領域分割部３３４は、ステップＳ１１０１で算出された対象組織抽出画像Ｍに基づいて、総位相分布Ｐ_{ｔｏｔａｌ}を、少なくとも２つ以上の領域、即ち、対象組織の所定領域と、対象組織の所定領域以外の領域と、に分割する。対象組織の所定領域とは、例えば、脳内部の画素値を１、脳外部の画素値を０とした画像において、画素値が１である領域と画素値が０である領域との境界部分である脳表領域である。また、対象組織の所定領域以外の領域とは、例えば、脳表領域を除いた脳内部の領域である。なお、脳表領域として定義される領域は、後述のステップＳ１１０６において、全領域局所位相分布算出部３３６が、対象組織の全領域における局所位相分布を算出する際に用いる球状カーネルρ_ｒの半径ｒ（カーネルサイズ）によって決定される。

ステップＳ１１０５において、背景位相分布算出部３３５は、ステップＳ１１０３で算出された総位相分布Ｐ_{ｔｏｔａｌ}から、対象組織の所定領域（例えば、脳表領域）における背景位相分布Ｐｅｄｇｅを算出する。以下、対象組織の所定領域における背景位相分布Ｐｅｄｇｅの算出法について説明する。

背景位相分布算出部３３５は、局所多項式近似処理により、総位相分布Ｐ_{ｔｏｔａｌ}から、対象組織の所定領域における背景位相分布Ｐｅｄｇｅを算出する。
背景位相分布算出部３３５は、局所領域内の任意画素の位置を（ｘ、ｙ、ｚ）とし、式（６）を満たす多項式係数Ａ_{ｋ、ｌ、ｎ}を算出する。

Ｋはｘ方向の多項式における次数、Ｌはｙ方向の多項式における次数、Ｎはｚ方向の多項式における次数をそれぞれ表す。
局所領域内の中心画素の位置を（ｘ_０、ｙ_０、ｚ_０）とすると、対象組織の所定領域における背景位相分布Ｐｅｄｇｅは、式（７）により算出できる。

背景位相分布算出部３３５は、対象組織の所定領域における画素の全てに対して、式（７）を実施することにより、背景位相分布Ｐｅｄｇｅを算出する。

ステップＳ１１０６において、全領域局所位相分布算出部３３６は、総位相分布Ｐ_{ｔｏｔａｌ}と、対象組織の所定領域における背景位相分布Ｐｅｄｇｅとを組み合わせて、対象組織の全領域における局所位相分布Ｐ’_{ｌｏｃａｌ}を算出する。以下、対象組織の全領域における局所位相分布Ｐ’_{ｌｏｃａｌ}の算出法について説明する。

全領域局所位相分布算出部３３６は、制約付き最小二乗処理により、対象組織の全領域における局所位相分布の推定値Ｐ’_{ｌｏｃａｌ}を算出する。

Ｗ_１は脳表領域における重みが小さい重み画像、Ｗ_２は脳表領域における重みが大きい重み画像をそれぞれ表す。
式（８）における第一項は、ＳＨＡＲＰ法と同様、脳内部の領域における、平均値の性質に基づく制約項を表している。平均値の性質における条件を満たす制約項に対しては、脳内部の領域における重みを大きくし、脳表領域における重みを小さくする。
式（８）における第二項は、脳表領域における背景位相分布に基づく制約項を表している。脳表領域における背景位相分布に基づく制約項に対しては、脳表領域における重みを大きくし、脳内部の領域における重みを小さくする。

ここで、図７を参照して、重み画像Ｗ_１及び重み画像Ｗ_２と式（８）との関係について説明する。図７（ａ）は、重み画像Ｗ_１の一次元プロファイルを示している。図７（ｂ）は、重み画像Ｗ_２の一次元プロファイルを示している。横軸は位置を示し、縦軸は重みを示している。

図７において、領域７０１と領域７０２とを合わせた領域は、脳表領域を示し、領域７０３は、脳表領域を除いた脳内部の領域（脳表領域以外の領域）を示している。同様に、領域７０４と領域７０５とを合わせた領域は、脳表領域を示し、領域７０６は、脳表領域を除いた脳内部の領域（脳表領域以外の領域）を示している。なお、領域７０１（領域７０４）及び領域７０２（領域７０５）は、球状カーネルρ_ｒの半径ｒよりも大きくなるように設定される。

図７に示すように、全領域局所位相分布算出部３３６は、領域７０１での重みが０となるように設定し、領域７０３での重みが１となるように設定し、領域７０１から領域７０３へと遷移する領域７０２での重みが０から１の間の値（線形に変化する値）となるように設定する。また、全領域局所位相分布算出部３３６は、領域７０４での重みが１となるように設定し、領域７０６での重みが０となるように設定し、領域７０４から領域７０６へと遷移する領域７０５での重みが０から１の間の値（線形に変化する値）となるように設定する。

即ち、全領域局所位相分布算出部３３６は、重み画像Ｗ_１及び重み画像Ｗ_２を用いて、脳内部の領域における重みを大きくし、脳表領域における重みを小さくすることで、式（８）に示す平均値の性質に基づく制約項を、ＳＨＡＲＰ法と同等の精度で算出することができる（式（５）参照）。
また、全領域局所位相分布算出部３３６は、重み画像Ｗ_１及び重み画像Ｗ_２を用いて、脳表領域における重みを大きくし、脳内部の領域における重みを小さくすることで、式（８）に示す脳表領域における背景位相分布に基づく制約項を、局所多項式近似処理により算出された脳表領域における背景位相分布に近づけることができる。
更に、全領域局所位相分布算出部３３６は、重み画像Ｗ_１及び重み画像Ｗ_２を用いて、領域７０２及び領域７０６における重みを、０から１の間で線形に変化させることで、脳表領域（領域７０１と領域７０２とを合わせた領域、領域７０４と領域７０５とを合わせた領域）と、脳内部の領域（領域７０３、領域７０６）とを、滑らかに接続することができる。

なお、全領域局所位相分布算出部３３６は、式（８）に示すような制約付き最小二乗処理のみならず、式（９）に示すような正則化付最小二乗処理（局所位相分布の正則化項を最小化する条件）を用いて、対象組織の全領域（脳表領域を含めた脳全体）における局所位相分布を算出してもよい。

例えば、全領域局所位相分布算出部３３６は、局所位相分布Ｐ_{ｌｏｃａｌ}の画素値のばらつきを考慮した正則化項を挿入して、式（９）を実施することにより、対象組織の全領域における局所位相分布を算出することも可能である。

λ_１は、対象組織の所定領域における制約の大きさを調整するための正則化パラメータである。λ_２は、局所位相分布Ｐ_{ｌｏｃａｌ}のばらつきを、ある程度許容するための正則化パラメータである。
式（９）における第一項は、ＳＨＡＲＰ法と同様、脳内部の領域における、平均値の性質に基づく制約項を表している。
式（９）における第二項は、脳表領域における背景位相分布に基づく制約項を表している。
式（９）における第三項は、位相分布から局所位相変化成分を過剰に除去してしまうことを防ぐ役割を持つ正則化項を表している。

式（８）及び式（９）より、本実施形態に係る計算機１０９の演算処理によれば、脳内部の領域における局所位相分布の算出精度は、従来のＳＨＡＲＰ法と略同様の算出精度を維持しつつ、脳表領域における局所位相分布の算出精度を、従来と比較して格段に向上させることができることが示された。即ち、本実施形態に係る計算機１０９の演算処理によれば、脳内部の領域における平均値の性質に基づく制約項及び脳表領域における背景位相分布に基づく制約項から、脳表領域を含めた脳全体の局所位相分布を高精度に算出できる
ことが示された。

［所望分布算出：Ｓ１００４］
Ｓ１００４において、計算機１０９は所望分布算出処理を行う。

所望分布算出部３４０は、ステップＳ１００３で算出された局所位相分布Ｐ_{ｌｏｃａｌ}を用いて、ＱＳＭ法やＳＷＩ法などにより、定量的磁化率分布、磁化率強調画像などの所望分布を算出する。例えば、所望分布算出部３４０は、局所磁場算出部と磁化率分布算出部とを備え、ＱＳＭ法により、定量的磁化率分布を算出する。また、例えば、所望分布算出部３４０は、マスク作成部と磁化率強調画像作成部とを備え、ＳＷＩ法により、磁化率強調画像を算出する。以下、所望分布算出部３４０が行う具体的な処理例として、ＱＳＭ法を用いた定量的磁化率分布の算出法と、ＳＷＩ法を用いた磁化率強調画像の算出法について説明する。

≪ＱＳＭ法を用いた定量的磁化率分布の算出法≫
まず、ＱＳＭ法による定量的磁化率分布の算出法について説明する。ＱＳＭ法は、局所位相分布Ｐ_{ｌｏｃａｌ}が、生体組織間の磁化率差を反映することを利用して、局所位相分布Ｐ_{ｌｏｃａｌ}から、局所的な磁場分布を算出し、磁場と磁化率の関係式を用いて定量的磁化率分布を算出する手法である。

位置ベクトルをｒとすると、生体組織間の磁化率差によって生じる相対的な磁場分布は、式（１０）で表される。

gはプロトンの核磁気回転比、Ｂ０は静磁場強度をそれぞれ表す。
ここで、局所的な磁場分布は、静磁場に関するマクスウェル方程式より、式（１１）で表される。

ここで、d（ｒ）は局所磁場分布、c（ｒ）は生体内の磁化率分布、aはベクトル（ｒ’−ｒ）と静磁場方向との為す角度、ｄ（ｒ）は点ダイポール磁場をそれぞれ表す。
即ち、式（１１）に示すように、局所磁場分布d（ｒ）は、磁化率分布c（ｒ）と点ダイポール磁場ｄ（ｒ）との畳み込み積分で表すことができる。従って、式（１１）の両辺をフーリエ変換することにより、式（１１）は、式（１２）に変換できる。

ｋ＝（ｋ_ｘ、ｋ_ｙ、ｋ_ｚ）は、ｋ空間上の位置ベクトルを表す。
D（ｋ）は、局所磁場分布d（ｒ）のフーリエ成分を表す。
Ｄ（ｋ）は、点ダイポール磁場ｄ（ｒ）のフーリエ成分を表す。
Ｘ（ｋ）は、生体内の磁化率分布c（ｒ）のフーリエ成分を表す。
式（１２）に示すように、生体内の磁化率分布c（ｒ）のフーリエ成分Ｘ（ｋ）は、局所磁場分布d（ｒ）のフーリエ成分D（ｋ）を、点ダイポール磁場ｄ（ｒ）のフーリエ成分Ｄ（ｋ）で除算することによって算出できる。
しかしながら、式（１２）は、Ｄ（ｋ）＝０近傍の領域において、その逆数が発散してしまうため、直接的にＸ（ｋ）を算出することができない。このＤ（ｋ）＝０近傍の領域はマジックアングル領域と呼ばれ、静磁場方向に対しておよそ５４．７度の２倍の頂角を持つ逆双円錐領域となる。

ＱＳＭ法は、マジックアングル領域の存在により不良条件逆問題（ill-conditioned inverse problem）に帰着されるため、不良条件逆問題を解くための幾つかの解法が提案されている。不良条件逆問題を解く代表的な方法としては、例えば、磁場と磁化率との関係式に基づく制約条件下で、磁場分布から算出した磁化率分布に対して平滑化処理を行うことを繰り返す方法、磁場分布及び点ダイポール磁場のｋ空間上の演算により磁化率分布を算出するＴＫＤ（Truncated-based K-space Division）法、ＴＫＤ法で算出した磁化率分布と閾値処理により微細構造を抽出した磁化率分布とを繰り返し演算により合成するIterative ＳＷＩＭ（Susceptibility Weighted Imaging and Mapping）法、正則化付最小二乗法を用いたＭＥＤＩ（Morphology enabled dipole inversion）法、などが挙げられる。なお、所望分布算出部３４０は、これらの何れの方法を用いてもよい。

≪ＳＷＩ法を用いた磁化率強調画像の算出法≫
次に、ＳＷＩ法による磁化率強調画像の算出法について説明する。ＳＷＩ法は、局所位相分布Ｐ_{ｌｏｃａｌ}から、磁化率分布を強調する磁化率強調マスクを作成し、計測した強度画像（絶対値画像）に、磁化率強調マスクを乗算することで、磁化率強調画像を算出する手法である。

所望分布算出部３４０は、所定の磁化率強調マスクを作成し、この磁化率強調マスクを、所定回数乗算する。その後、所望分布算出部３４０は、任意のエコー時間の絶対値画像に、当該磁化率強調マスクを乗算することで、磁化率強調画像を算出する。磁化率強調画像とは、磁化率の高い所が暗くなり、磁化率の低い所が明るくなる画像である。操作者は、当該画像に基づいて、磁化率が高くなっているか、或いは低くなっているかを判断することができる。なお、所望分布算出部３４０は、公知のさまざまな手法を用いて、磁化率強調画像を算出することが可能であり、ＳＷＩ法に限定されるものではない。

［画像表示：Ｓ１００５］
Ｓ１００５において、計算機１０９は画像表示処理を行う。

計算機１０９は、ステップＳ１００４で算出された定量的磁化率分布や磁化率強調画像を、表示部１１０のディスプレイ画面上に表示させる。また、計算機１０９は、ステップＳ１００４で算出された各種データを外部記憶装置１１１或いは記憶装置に記憶させ、所望の表示装置に、定量的磁化率分布や磁化率強調画像を表示させることも可能である。

ここで、図８を参照して、ＭＲＩ装置１００の表示部１１０に表示される画像の一例について説明する。

図８（ａ）は、本実施形態に係る計算機１０９によって演算処理が施された画像（磁化率強調画像）であり、図８（ｂ）は、図８（ａ）の領域８０１における拡大図である。図８（ｃ）は、従来の計算機によって演算処理が施された画像（磁化率強調画像）であり、図８（ｄ）は、図８（ｃ）の領域８０２における拡大図である。

図８（ａ）及び図８（ｂ）に示す画像と図８（ｃ）及び図８（ｄ）に示す画像とを比較すると、図８（ａ）及び図８（ｂ）に示す画像では、脳表領域が鮮明に撮像されている一方で、図８（ｃ）及び図８（ｄ）に示す画像では、脳表領域が太い実線で囲まれたように撮像されていることがわかる。即ち、本実施形態に係る計算機１０９は、従来の計算機と比較して、太い実線で囲まれた領域の分だけ、より多くの情報を算出できることがわかる。

従って、操作者は、脳全体における生体組織の構造や組成などを診断する場合、図８（ｃ）及び図８（ｄ）に示す画像を用いるよりも、図８（ａ）及び図８（ｂ）に示す画像を用いることで、より適切な診断を行うことができる。更に、操作者は、くも膜下出血や脳梗塞など、脳の表面に発生した病気を診断する場合、図８（ａ）及び図８（ｂ）に示す画像を用いることで、正確な診断を行うことができる。即ち、本実施形態に係る計算機１０９が所定の演算処理を行うことで、操作者は、脳表領域において信頼性の高い画像を取得することができる。

以上、上述の各処理を行うことで、本実施形態に係る計算機１０９が行う演算処理が終了する。なお、上述の図４及び図６に示すフローは一例であり、その一部を省略することや、別の処理を追加することが可能であることは勿論である。

本実施形態に係るＭＲＩ装置によれば、計算機１０９が、複素画像から、例えば脳領域を抽出し、総位相分布を脳表領域と脳内部の領域とに分割して、脳表領域における背景位相分布を算出し、総位相分布と該背景位相分布とを組み合わせることで、脳全体の局所位相分布を算出する。これにより、脳表領域を含めた脳全体の局所位相分布を高精度に算出するＭＲＩ装置を実現することができる。

更に、本実施形態に係るＭＲＩ装置によれば、高精度な局所位相分布を用いて、定量的磁化率分布や磁化率強調画像を算出することができる。これにより、操作者は、鮮明な領域が拡張された画像を取得することができるため、ＭＲＩ装置の信頼性を高められる。

１００，１２０，１３０ 磁気共鳴イメージング装置
１０１ 被検体
１０２ 静磁場生成磁石
１０３ 傾斜磁場印加部
１０５ 送信部
１０６ 受信部
１０９ 計算機
１１０ 表示部
３３１ 対象組織抽出部
３３２ 位相分布算出部
３３３ 位相折り返し補正部
３３４ 領域分割部
３３５ 背景位相分布算出部
３３６ 全領域局所位相分布算出部（局所位相分布算出部）

本発明は、磁気共鳴イメージング（MRI:magnetic resonance imaging）技術に関する。特に、Gradient Echo法で撮像された複素画像を用いて、生体組織間の磁化率差に起因した局所位相分布を算出する画像処理技術に関する。

磁気共鳴イメージング装置は、静磁場内に配置される被検体の水素原子核（プロトン）が、特定周波数の高周波磁場に共鳴する核磁気共鳴（NMR:nuclear magnetic resonance）現象を利用した非侵襲な医用画像診断装置である。被検体から発生する核磁気共鳴信号は、プロトン密度や緩和時間など様々な物性値によって変化するため、磁気共鳴イメージング装置は、核磁気共鳴信号に基づいて、生体組織の構造や組成、細胞性状などの様々な生体情報を描出できる。

近年、磁気共鳴イメージング装置により測定可能な物性値として、生体組織間の磁化率差が注目されている。磁化率とは、静磁場中の物質が磁気分極（磁化）する度合いを表した物性値である。生体内には、静脈血中のデオキシヘモグロビンや鉄たんぱく質などの常磁性体と、生体組織の大部分を占める水や石灰化の基となるカルシウムなどの反磁性体とが存在する。磁気共鳴イメージング装置は、生体組織間の磁化率差を画像化することができるため、脳虚血疾患の診断やがんの放射線治療効果の予測、神経変性疾患の鑑別、などへの可能性が期待されている。生体組織間の磁化率差を画像化する方法としては、例えば、定量的磁化率マッピング（QSM:quantitative susceptibility mapping）法、磁化率強調画像化（SWI:Susceptibility Weighted Imaging）法、などが知られている。

磁気共鳴イメージング装置は、生体組織間の磁化率差に起因する局所位相分布（局所的な位相変化の分布）を算出し、当該局所位相分布を用いてＱＳＭ法により定量的磁化率分布を算出する、或いは、当該局所位相分布を用いてＳＷＩ法により磁化率強調画像を算出する。磁気共鳴イメージング装置は、当該局所位相分布を算出する際、位相折り返し補正処理を行って総位相分布を算出し、総位相分布から被検体の形状などに起因して生じる背景位相分布を除去する背景位相除去処理を行う。背景位相除去処理における代表的な方法として、以下の二つの方法が提案されている。

一つ目の方法は、ＰＤＦ（Projection onto Dipole Fieled）法と呼ばれる方法である（例えば、特許文献１参照）。ＰＤＦ法は、計測した脳内部の位相分布から脳外部の磁化率分布を推定し、該磁化率分布によって生じる脳内部の位相分布を背景位相分布として算出し、総位相分布から背景位相分布を除去することで局所位相分布を算出する。

二つ目の方法は、ＳＨＡＲＰ（Sophisticated Harmonic Artifact Reduction for Phase date）法と呼ばれる方法である（例えば、非特許文献１参照）。ＳＨＡＲＰ法は、背景位相分布が球面調和関数近似できることを利用する。即ち、任意の点における値が、その点を中心とした任意の半径を有する球（カーネル）内の点における値の平均に等しい、という平均値の性質（Spherical mean value property）を用いて、総位相分布から背景位相分布と局所位相分布とを分離する。

米国特許第９４４８２８９号明細書

Ｆｅｒｄｉｎａｎｄ Ｓｃｈｗｅｓｅｒ他、 "Ｑｕａｎｔｉｔａｔｉｖｅ ｉｍａｇｉｎｇ ｏｆ ｉｎｔｒｉｎｓｉｃ ｍａｇｎｅｔｉｃ ｔｉｓｓｕｅ ｐｒｏｐｅｒｔｉｅｓ ｕｓｉｎｇ ＭＲＩ ｓｉｇｎａｌ ｐｈａｓｅ： Ａｎ ａｐｐｒｏａｃｈ ｔｏ ｉｎ ｖｉｖｏ ｂｒａｉｎ ｉｒｏｎ ｍｅｔａｂｏｌｉｓｍ？"、 Ｎｅｕｒｏｉｍａｇｅ、２０１１年、５４巻、２７８９−２８０７頁

しかしながら、特許文献１に記載の方法では、脳外部の推定磁化率で再現できない部分（例えば、副鼻腔など）の背景位相分布を除去することができないという問題がある。また、非特許文献１に記載の方法では、脳内部と脳外部との境界部分（脳表領域）における局所位相分布の算出精度が低下してしまうという問題がある。

本発明は、上記した課題を解決するためになされたものであり、脳表領域を含めた脳全体の局所位相分布を高精度に算出する画像処理方法、画像処理プログラム及び計算機を提供することを課題とする。

前記課題を解決するために、本開示の実施形態に係る計算機は、磁気共鳴イメージングで撮像した被検体の複素画像から前記被検体の対象組織を抽出した対象組織抽出画像を算出する対象組織抽出部と、前記複素画像から位相分布を算出する位相分布算出部と、前記位相分布の位相折り返しを補正する位相折り返し補正部と、前記対象組織抽出画像に基づいて、補正後の位相分布を、前記対象組織抽出画像における対象組織の表面領域と、前記表面領域以外の前記対象組織の内部領域とに分割する領域分割部と、前記表面領域の背景位相分布を算出する背景位相分布算出部と、前記補正後の位相分布と前記背景位相分布とを組み合わせて、前記表面領域における局所位相分布を算出する局所位相分布算出部と、を備えたことを特徴とする。

また本開示の実施形態に係る画像処理方法は、磁気共鳴イメージングで撮像した被検体の複素画像を用いて、前記被検体の対象組織の所定領域における局所位相分布を算出する画像処理方法であって、上記計算機の各部により実現される処理を含むものであり、例えば、核磁気共鳴信号によって構成された前記複素画像から前記被検体の対象組織を抽出した対象組織抽出画像を算出するステップと、前記複素画像から位相分布を算出するステップと、前記位相分布の位相折り返しを補正するステップと、前記対象組織抽出画像に基づいて、補正後の位相分布を、前記対象組織の表面領域と前記所定領域以外の前記対象組織の内部領域とに分割するステップと、前記表面領域の背景位相分布を算出するステップと、を含み、前記補正後の位相分布と前記背景位相分布とを組み合わせて、前記表面領域における局所位相分布を算出することを特徴とする。

さらに本開示の実施形態の画像処理プログラムは、上記画像処理方法を計算機に機能させるプログラムである。

本発明によれば、脳表領域を含めた脳全体の局所位相分布を高精度に算出することができる。

（ａ）は、本実施形態に係る垂直磁場方式の磁気共鳴イメージング装置の外観図、（ｂ）は、本実施形態に係る水平磁場方式の磁気共鳴イメージング装置の外観図、（ｃ）は、本実施形態に係る開放感を高めた磁気共鳴イメージング装置の外観図である。 本実施形態に係る磁気共鳴イメージング装置の概略構成を示すブロック図である。 本実施形態に係る計算機の機能ブロック図である。 本実施形態に係る計算機における演算処理の流れを説明するフローチャートである。 ＲＳＳＧ（RF-spoiled-Steady-state Acquisition with Rewound Gradient-Echo）シーケンスのパルスシーケンスを説明するための図である。 本実施形態に係る局所位相分布算出部における局所位相分布算出処理の流れを説明するフローチャートである。 （ａ）は、本実施形態に係る全領域局所位相分布算出部で用いる重み画像Ｗ_１を説明するための図、（ｂ）は、本実施形態に係る全領域局所位相分布算出部で用いる重み画像Ｗ_２を説明するための図である。 （ａ）は、本実施形態に係る計算機によって演算処理が施された画像を示す図、（ｂ）は、（ａ）の一部領域における拡大図、（ｃ）は、従来の計算機によって演算処理が施された画像を示す図、（ｄ）は、（ｃ）の一部領域における拡大図である。

本発明を適用する実施形態の一例を説明する。以下、実施形態を説明するための全図において、特に断らない限り、同一機能を有するものは同一符号を付し、繰り返しの説明は省略する。また、以下の記述により本発明が限定されるものではない。

≪ＭＲＩ装置の全体構成≫
図１は、本実施形態に係るＭＲＩ装置の外観図である。図１（ａ）は、ソレノイドコイルで静磁場を生成するトンネル型磁石を用いた水平磁場方式のＭＲＩ装置１００である。図１（ｂ）は、開放感を高めるために磁石を上下に分離したオープン型の垂直磁場方式のＭＲＩ装置１２０である。また、図１（ｃ）は、図１（ａ）と同じトンネル型磁石を用い、磁石の奥行を短くし、かつ、斜めに傾けることによって、開放感を高めたＭＲＩ装置１３０である。なお、本実施形態において、ＭＲＩ装置の構成は、特に限定されるものではなく、図１に示す何れの構成であってもよいし、公知の構成であってもよい。

以下、特に区別する必要がない場合は、ＭＲＩ装置１００で代表する。ＭＲＩ装置１００において、ＭＲＩ装置１００の静磁場方向をｚ方向、ｚ方向に垂直な２方向のうち、測定対称である被検体１０１を載置するベッド面に平行な方向をｘ方向、他方向をｙ方向とする。

図２は、本実施形態に係るＭＲＩ装置の概略構成を示すブロック図である。図２に示すように、ＭＲＩ装置１００は、被検体１０１が配置される空間に静磁場を生成する静磁場生成磁石１０２と、被検体１０１に高周波磁場を送信する送信用高周波コイル（送信部）１０５と、高周波磁場の送信により被検体１０１から発生する核磁気共鳴信号を受信する受信用高周波コイル（受信部）１０６と、核磁気共鳴信号に位置情報を付加する傾斜磁場を印加する傾斜磁場印加部１０３と、送信部１０５、受信部１０６、及び傾斜磁場印加部１０３を制御するとともに、核磁気共鳴信号に対して演算処理を行う計算機１０９と、計算機１０９により演算処理された画像を表示する表示部１１０と、を備える。

更に、ＭＲＩ装置１００は、シムコイル１０４、送信機１０７、受信機１０８、外部記憶装置１１１、傾斜磁場用電源部１１２、シム用電源部１１３、シーケンス制御装置１１４、入力装置１１５、などを備える。

静磁場生成磁石１０２は、被検体１０１が配置される空間に静磁場を生成する磁石であり、ＭＲＩ装置の構成に応じて、種々の形態を採用することができる。傾斜磁場印加部１０３は、核磁気共鳴信号に空間的な位置情報を付加するために、被検体１０１が配置される空間において、ｘ方向、ｙ方向、ｚ方向の各方向に傾斜磁場を印加するコイルである。各方向に印加される傾斜磁場（ｘ方向傾斜磁場、ｙ方向傾斜磁場、ｚ方向傾斜磁場）は、例えば、位相エンコード傾斜磁場、リードアウト傾斜磁場、スライス傾斜磁場である。

シムコイル１０４は、被検体１０１が配置される空間に生成される静磁場の不均一性を打ち消して除去するなど、静磁場分布を調整するために使用される電流コイルである。送信部１０５は、被検体１０１の水素原子核を励起させる高周波磁場を、被検体１０１に送信する送信用高周波コイルである。受信部１０６は、高周波磁場の送信により被検体１０１から発生する核磁気共鳴信号を受信する受信用高周波コイルである。なお、送信部１０５と受信部１０６とは、別々の構成としてもよいし、兼用する構成としてもよい。

送信機１０７は、シーケンス制御装置１１４によって制御され、送信部１０５へと送信する高周波磁場を生成する。受信機１０８は、シーケンス制御装置１１４によって制御され、受信部１０６によって受信された核磁気共鳴信号を検波し、複素信号を計算機１０９へと送信する。なお、送信部１０５と受信部１０６とは、別々の構成としてもよいし、兼用する構成としてもよい。

計算機１０９は、ＣＰＵ（Central Processing Unit）、メモリ、記憶装置などを備える情報処理装置であり、計算機１０９には、表示部１１０、外部記憶装置１１１、入力装置１１５などが接続される。計算機１０９は、記憶装置に記憶される制御プログラムを読み出して、ワークエリアに展開し、当該制御プログラムを実行することで、ＭＲＩ装置１００全体の動作を制御するとともに、各種演算処理を行う。なお、計算機１０９は、ＭＲＩ装置１００の一部であってもよいし、ＭＲＩ装置１００から独立した情報処理装置であって、ＭＲＩ装置１００とデータの送受信が可能な情報処理装置であってもよい。更に、計算機１０９が実現する各種機能のうち、全部または一部の機能は、ＡＳＩＣ（Application Specific Integrated Circuit）、ＦＰＧＡ（field-programmable gate array）などのハードウェアによって実現されていてもよい。

表示部１１０は、例えば、液晶ディスプレイ、有機ＥＬディスプレイ、等で構成される。表示部１１０は、計算機１０９から取得したデータに基づいて、ディスプレイ画面上に、演算結果、各種画像などを表示する。表示部１１０は、例えば、後述する局所位相分布、定量的磁化率分布、磁化率強調画像などをディスプレイ画面上に表示する。

外部記憶装置１１１は、記憶装置とともに、計算機１０９が行う各種演算処理に使用されるデータ、演算処理中に生成されるデータ、各種演算処理により取得されるデータ、入力装置１１５を介して入力される各種条件、各種パラメータ、などを記憶する。

傾斜磁場用電源部１１２は、傾斜磁場印加部１０３の駆動用電源であり、シーケンス制御装置１１４からの制御指示に基づいて、傾斜磁場印加部１０３の駆動を制御する。シム用電源部１１３は、シムコイル１０４の駆動用電源であり、シーケンス制御装置１１４からの制御指示に基づいて、シムコイル１０４の駆動を制御する。

シーケンス制御装置１１４は、送信機１０７、受信機１０８、傾斜磁場用電源部１１２、シム用電源部１１３、などＭＲＩ装置１００の各部の動作を制御する。例えば、シーケンス制御装置１１４は、外部記憶装置１１１や記憶装置に記憶される設定情報（例えば、後述のパルスシーケンスと呼ばれるタイムチャート）に基づいて、高周波磁場の印加タイミング、傾斜磁場（ｘ方向傾斜磁場、ｙ方向傾斜磁場、ｚ方向傾斜磁場）の印加タイミング、核磁気共鳴信号の受信タイミング、などを制御する。

入力装置１１５は、操作者による入力操作を受け付ける機能を備えており、例えば、キーボード、カメラ、マウス、タッチパネル、音声入力を受け付けるマイク、などで構成される。入力装置１１５は、計算機１０９が行う各種演算処理や計測処理などに必要な各種条件、各種パラメータ、などを操作者が入力するためのインタフェースである。なお、各種条件、各種パラメータ、などは、予め設定しておくことも可能であるし、操作者が入力装置１１５を介して適宜設定することも可能である。

≪計算機の構成≫
次に、図３を参照して、本実施形態に係る計算機１０９の構成について説明する。本実施形態では対象組織として被検体１０１の脳を適用する場合を一例に挙げて説明する。

ここで、計算機１０９の構成の説明に先立ち、ＳＨＡＲＰ法を用いた一般的な背景位相除去処理の概要について説明する。

まず、ＭＲＩ装置は、Gradient Echo法を用いて、任意のエコー時間（高周波磁場を送信してから核磁気共鳴信号を受信するまでの時間）で撮像された複素画像から位相分布Ｐを算出する。そして、ＭＲＩ装置は、位相分布Ｐに生じている−π〜+πの範囲に折り返された位相を補正する位相折り返し補正処理を行う。

位相折り返し補正処理後の位相分布（総位相分布）Ｐ_{ｔｏｔａｌ}は、対象組織の形状などに起因する背景位相分布Ｐ_ｂｋｇｒ、生体組織間の磁化率差に起因する局所位相分布Ｐ_{ｌｏｃａｌ}を用いて、式（１）で表される。

そして、ＭＲＩ装置は、任意画素における画素値と、その画素を中心とした所定の半径を有する球(カーネル)内の画素における画素値の平均とが略等しいという平均値の性質を用いて、総位相分布Ｐ_{ｔｏｔａｌ}から、背景位相分布Ｐ_ｂｋｇｒと局所位相分布Ｐ_{ｌｏｃａｌ}とを分離する。

背景位相分布Ｐ_ｂｋｇｒは、平均値の性質によって、式（２）を満たす。

δはデルタ関数、ρ_ｒは半径ｒの球状カーネル、符号（〇で囲った×）は畳み込み積分をそれぞれ表す。また、Ｍ_{ｓｈｒｉｎｋ}は、被検体１０１の対象組織を抽出した対象組織抽出画像Ｍから、球状カーネルρ_ｒの半径ｒに応じて縮小されたマスクを表す。

式（２）で表される畳み込み積分は、フーリエ変換行列Ｆ、逆フーリエ変換行列Ｆ^−１を用いることで、式（３）に変換できる。

但し、Ｃ＝Ｆ(δ−ρ_ｒ)

式（３）に式（１）を代入することで、式（４）の関係式が得られる。

そして、ＭＲＩ装置は、総位相分布Ｐ_{ｔｏｔａｌ}から、式（４）の関係式を満たす局所位相分布Ｐ_{ｌｏｃａｌ}を算出し、最小二乗処理によって、式（５）で表される局所位相分布Ｐ’_{ｌｏｃａｌ}を算出する。

即ち、ＳＨＡＲＰ法を用いた一般的な背景位相除去処理では、マスクＭ_{ｓｈｒｉｎｋ}の範囲における局所位相分布しか算出することができず、脳内部と脳外部との境界部分（脳表領域）における局所位相分布の算出精度が低下するという問題が生じていた。

そこで、本実施形態に係るＭＲＩ装置１００（計算機１０９）は、総位相分布を、脳表領域と脳内部の領域とに領域分割し、総位相分布と脳表領域における背景位相分布とを組み合わせることで、脳表領域を含めた脳全体の局所位相分布を算出する。これにより、脳表領域における局所位相分布の算出精度を向上させることができる。以下、これを実現する計算機１０９の機能構成について説明する。

計算機１０９は、計測制御部３１０と、画像再構成部３２０と、局所位相分布算出部３３０と、所望分布算出部３４０と、を備える。更に、局所位相分布算出部３３０は、対象組織抽出部３３１、位相分布算出部３３２、位相折り返し補正部３３３、領域分割部３３４、背景位相分布算出部３３５、全領域局所位相分布算出部３３６と、を備える。

計測制御部３１０は、例えば、Gradient Echo法を用いて、送信部１０５から被検体１０１への高周波磁場の送信により、被検体１０１から発生する任意のエコー時間の核磁気共鳴信号を複素信号として計測する。計測制御部３１０は、予め設定されたパルスシーケンスに従って、シーケンス制御装置１１４を制御することで、高周波磁場の印加タイミング、傾斜磁場の印加タイミング、核磁気共鳴信号の受信タイミング、などを制御する。なお、計測制御部３１０は、血流などの流れの影響を補償する流速補正（Flow Compensation）傾斜磁場パルスを、ｘ方向、ｙ方向、ｚ方向の各方向に印加してもよい。

画像再構成部３２０は、計測制御部３１０によって計測された複素信号（核磁気共鳴信号）を、例えば、ｋｒ軸、ｋｐ軸、ｋｓ軸を座標軸とする３次元のｋ空間に配置してフーリエ変換することにより、任意画素の画素値が複素数で表される複素画像を再構成する。画像再構成部３２０は、ｋ空間に核磁気共鳴信号を充填する際、ｋ空間の座標軸（例えば、ｋｒ軸）に平行に、１行ずつ核磁気共鳴信号を充填していく方法を用いてもよいし（Cartesian充填法）、ｋ空間の原点を中心に、角度を変化させながら放射状に核磁気共鳴信号を充填していく方法を用いてもよい（Radial Scan充填法）。

局所位相分布算出部３３０は、画像再構成部３２０によって再構成された複素画像に基づいて、対象組織の全領域における局所位相分布を算出する。具体的には、対象組織抽出部３３１が、複素画像から被検体１０１の対象組織を抽出した対象組織抽出画像Ｍを生成する。そして、位相分布算出部３３２が、複素画像から位相分布Ｐを算出する。そして、位相折り返し補正部３３３が、位相分布Ｐの位相折り返しを補正し、総位相分布Ｐ_{ｔｏｔａｌ}を算出する。そして、領域分割部３３４が、対象組織抽出画像Ｍに基づいて、総位相分布Ｐ_{ｔｏｔａｌ}を少なくとも２つ以上の領域（例えば、所定領域、所定領域以外の領域）に分割する。そして、背景位相分布算出部３３５が、所定領域における背景位相分布Ｐｅｄｇｅを算出する。そして、全領域局所位相分布算出部３３６が、総位相分布Ｐ_{ｔｏｔａｌ}と、所定領域における背景位相分布Ｐ_ｅｄｇｅとを組み合わせることで、所定領域を含めた対象組織の全領域における局所位相分布Ｐ’_{ｌｏｃａｌ}を算出する。

所望分布算出部３４０は、局所位相分布算出部３３０によって算出された局所位相分布Ｐ’_{ｌｏｃａｌ}を用いて、定量的磁化率分布や磁化率強調画像などの所望分布を算出し、当該所望分布を表示部１１０へと出力する。

上述のように、本実施形態に係る計算機１０９は、複素画像から対象組織のみを抽出し、総位相分布を対象組織の所定領域と所定領域以外の領域とに分割して、所定領域における背景位相分布を算出し、当該背景位相分布を利用して、対象組織の全領域における局所位相分布を算出する。これにより、カーネルの半径に応じて脳表領域で誤差が大きくなり、脳表領域における局所位相分布の算出精度が低下するという従来のＳＨＡＲＰ法に生じていた問題点を回避し、脳表領域を含めた脳全体における局所位相分布を高精度に算出することができる。

≪計算機における演算処理の流れ≫
次に、図４乃至図８を参照して、本実施形態に係る計算機１０９が行う演算処理の流れについて説明する。図４は、本実施形態に係る計算機１０９における演算処理の流れを説明するフローチャートである。図５は、ＲＳＳＧシーケンスのパルスシーケンスを説明するための図である。図６は、本実施形態に係る局所位相分布算出部３３０における局所位相分布算出処理の流れを説明するフローチャートである。図７は、本実施形態に係る全領域局所位相分布算出部３３６で用いる重み画像を説明するための図である。図８は、本実施形態に係る表示部１１０に表示される画像の一例を示す図である。

［計測：Ｓ１００１］
ステップＳ１００１において、計算機１０９は計測処理を行う。

操作者が、入力装置１１５を介して各種パラメータ（例えば、核磁気共鳴信号の数、エコー時間、エコー時間の間隔など）を設定し、計算機１０９に撮像開始指示が入力されると、計測制御部３１０は、予め設定されたパルスシーケンスに従って、シーケンス制御装置１１４を制御する。そして、計測制御部３１０は、任意のエコー時間の核磁気共鳴信号を複素信号として計測する。

ここで、図５を参照して、計測制御部３１０が計測に使用するGradient Echo型のパルスシーケンスについて説明する。図５に示すＲＳＳＧ（RF-spoiled-Steady-state Acquisition with Rewound Gradient-Echo）シーケンス５１０は、生体組織間の磁化率差によって生じる局所的な磁場変化に鋭敏なパルスシーケンスである。なお、計測制御部３１０が計測に使用するシーケンスは、ＲＳＳＧシーケンス５１０に限定されるものではない。

図５において、ＲＦは高周波磁場を示し、Ｇｓはスライス傾斜磁場を示し、Ｇｐは位相エンコード傾斜磁場を示し、Ｇｒはリードアウト傾斜磁場を示している。また、ＴＥはエコー時間を示し、ＴＲは繰り返し時間を示している。更に、５０１はスライス傾斜磁場パルス、５０２は高周波磁場パルス、５０３、５０８はスライスエンコード傾斜磁場パルス、５０４、５０９は位相エンコード傾斜磁場パルス、５０５、５０６はリードアウト傾斜磁場パルス、５０７は核磁気共鳴信号を示している。なお、ハイフン以下の数字は、計測制御部３１０が行う繰り返しの計測が何回目であるかを示しており、例えば、「−２」は、計測制御部３１０が行う繰り返しの計測が２回目であることを示している。

計測制御部３１０は、ＲＳＳＧシーケンス５１０に従って、シーケンス制御装置１１４を制御することで、スライス傾斜磁場パルス５０１、高周波磁場パルス５０２、スライスエンコード傾斜磁場パルス５０３、位相エンコード傾斜磁場パルス５０４、リードアウト傾斜磁場パルス５０５、リードアウト傾斜磁場パルス５０６、スライスエンコード傾斜磁場パルス５０８、位相エンコード傾斜磁場パルス５０９、などの印加タイミングを適切に制御することで、画像再構成に必要なデータを計測する。

図５に示すＲＳＳＧシーケンス５１０では、繰り返し時間ＴＲ内に、１つの核磁気共鳴信号５０７を計測する。
まず、静磁場内に被検体１０１が配置されると、静磁場強度に応じて被検体１０１内の水素原子核が、所定の周波数を持つため、送信部１０５は、被検体１０１の計測領域に、当該周波数と一致するような高周波磁場パルス５０２−１を送信し、被検体１０１の水素原子核を励起させる。この際、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｚ方向に、スライス傾斜磁場パルス５０１−１を印加し、被検体１０１における所定のスライス（３次元画像の撮像位置）での磁化を選択励起させる。

なお、計測制御部３１０は、繰り返しの計測を行う際、高周波磁場パルス５０２の位相を、所定角度（例えば、１１７度、１２２度）ずつ変更する。例えば、繰り返しの計測が１回目の場合、所定角度を１１７度とし、繰り返しの計測が２回目の場合、角度を１１７×２度とし、繰り返しの計測が３回目の場合、角度を１１７×３度とする。これにより、次の高周波磁場パルスを印加する際、余分な信号成分を残すことなく適切な計測を行うことが可能になる。

続いて、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｚ方向に、スライスエンコード傾斜磁場パルス５０３−１を印加して、磁化の位相に、スライスエンコード方向の位置情報を付加する。また、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｘ方向に、位相エンコード傾斜磁場パルス５０４−１を印加して、磁化の位相に、位相エンコード方向の位置情報を付加する。即ち、傾斜磁場印加部１０３は、核磁気共鳴信号の読み出しに先だって、ｚ方向及びｘ方向における磁化の位相に、位置情報を付加する。

スライスエンコード傾斜磁場パルス５０３−１は、繰り返し時間ＴＲ毎に、磁場強度が変化する。例えば、図５に示す７本の横線は、異なる強さの傾斜磁場を示しており、上方向の矢印に従って、正の傾斜磁場が強くなることを示している。同様に、位相エンコード傾斜磁場パルス５０４−１も、繰り返し時間ＴＲ毎に、磁場強度が変化する。例えば、図５に示す７本の横線は、異なる強さの傾斜磁場を示しており、上方向の矢印に従って、正の傾斜磁場が強くなることを示している。

なお、図５では、スライスエンコード傾斜磁場パルス５０３−１及び位相エンコード傾斜磁場パルス５０４−１における磁場強度が、７通りに変化する場合を一例に挙げて説明しているが、これに限定されるものではない。

続いて、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｙ方向に、ディフェーズ用のリードアウト傾斜磁場パルス５０５−１を印加して、画素内における磁化の位相を分散させる。そして、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｙ方向に、リードアウト傾斜磁場パルス５０６−１を印加して、磁化の位相に、リードアウト方向の位置情報を付加する。これにより、計測制御部３１０は、繰り返し時間ＴＲ内に、１つの核磁気共鳴信号５０７を計測することが可能になる。

最後に、傾斜磁場印加部１０３は、被検体１０１の計測領域におけるｚ方向に、リフェーズ用のスライスエンコード傾斜磁場パルス５０８−１を印加し、更に、被検体１０１の計測領域におけるｘ方向に、リフェーズ用の位相エンコード傾斜磁場パルス５０９−１を印加して、ディフェーズされた画素内における磁化の位相を収束させる。

スライスエンコード傾斜磁場パルス５０８−１は、繰り返し時間ＴＲ毎に、磁場強度が変化する。例えば、図５に示す７本の横線は、異なる強さの傾斜磁場を示しており、下方向の矢印に従って、負の傾斜磁場が強くなることを示している。同様に、位相エンコード傾斜磁場パルス５０９−１も、繰り返し時間ＴＲ毎に、磁場強度が変化する。例えば、図５に示す７本の横線は、異なる強さの傾斜磁場を示しており、下方向の矢印に従って、負の傾斜磁場が強くなることを示している。

計測制御部３１０は、上述の手順を、高周波磁場パルス５０２の位相、スライスエンコード傾斜磁場パルス５０３、５０８の強度、位相エンコード傾斜磁場パルス５０４、５０９の強度、などを変化させながら繰り返し行い、画像再構成に必要とされる核磁気共鳴信号５０７を計測する。なお、ＲＳＳＧシーケンス５１０により得られる複素画像の絶対値成分は、エコー時間ＴＥが短く設定されるとＴ１(縦緩和時間)強調画像となり、エコー時間ＴＥが長く設定されると画素内の位相分散を反映したＴ２^＊（見かけの横緩和時間）強調画像となる。本実施形態では、エコー時間ＴＥを長く設定し、Ｔ２^＊強調画像の撮像を行う。

［画像再構成：Ｓ１００２］
ステップＳ１００２において、計算機１０９は画像再構成処理を行う。

画像再構成部３２０は、ステップＳ１００１で計測された核磁気共鳴信号を、ｋ空間上に配置してフーリエ変換することにより、複素画像を再構成する。撮像方法としては、特に限定されるものではなく、カーテシアン撮像、ノンカーテシアン撮像、マルチエコーエコープラナーイメージング法、などの公知の方法を適用することができる。

［局所位相分布算出処理：Ｓ１００３］
ステップＳ１００３において、計算機１０９は局所位相分布算出処理を行う。

ここで、図６を参照して、局所位相分布算出処理の流れについて詳細に説明する。局所位相分布算出部３３０は、ステップＳ１００２で再構成された複素画像に基づいて、対象組織の全領域における局所位相分布を算出する。

ステップＳ１１０１において、対象組織抽出部３３１は、画像再構成部３２０によって構成された複素画像の絶対値成分（絶対値画像）における画素値の大きさに基づいて、被検体１０１の対象組織を抽出した対象組織抽出画像Ｍを算出する。対象組織抽出部３３１は、例えば、判別分析法などの公知の方法を用いて、絶対値画像の画素値のヒストグラムから、ノイズ成分と信号成分とを分離し、被検体１０１の対象組織を抽出した対象組織抽出画像Ｍ（例えば、脳内部の画素値を１、脳外部の画素値を０とした画像）を算出する。なお、対象組織抽出部３３１は、例えば、領域拡大（Region Growing）法などの公知の方法を用いて、複素画像から磁化率分布の算出に不要な領域（例えば、皮下脂肪などの領域）を除去して、磁化率分布の算出に必要な領域である対象組織の実質領域のみを抽出した対象組織抽出画像Ｍを算出してもよい。

ステップＳ１１０２において、位相分布算出部３３２は、画像再構成部３２０によって再構成された複素画像の偏角成分に基づいて、位相分布Ｐを算出する。

ステップＳ１１０３において、位相折り返し補正部３３３は、例えば、領域拡大法など公知の方法を用いて、ステップＳ１１０２で算出された位相分布Ｐに生じている−π〜＋πの範囲に折り返された位相を補正し、総位相分布Ｐ_{ｔｏｔａｌ}を算出する。

ステップＳ１１０４において、領域分割部３３４は、ステップＳ１１０１で算出された対象組織抽出画像Ｍに基づいて、総位相分布Ｐ_{ｔｏｔａｌ}を、少なくとも２つ以上の領域、即ち、対象組織の所定領域と、対象組織の所定領域以外の領域と、に分割する。対象組織の所定領域とは、例えば、脳内部の画素値を１、脳外部の画素値を０とした画像において、画素値が１である領域と画素値が０である領域との境界部分である脳表領域である。また、対象組織の所定領域以外の領域とは、例えば、脳表領域を除いた脳内部の領域である。なお、脳表領域として定義される領域は、後述のステップＳ１１０６において、全領域局所位相分布算出部３３６が、対象組織の全領域における局所位相分布を算出する際に用いる球状カーネルρ_ｒの半径ｒ（カーネルサイズ）によって決定される。

ステップＳ１１０５において、背景位相分布算出部３３５は、ステップＳ１１０３で算出された総位相分布Ｐ_{ｔｏｔａｌ}から、対象組織の所定領域（例えば、脳表領域）における背景位相分布Ｐｅｄｇｅを算出する。以下、対象組織の所定領域における背景位相分布Ｐｅｄｇｅの算出法について説明する。

背景位相分布算出部３３５は、局所多項式近似処理により、総位相分布Ｐ_{ｔｏｔａｌ}から、対象組織の所定領域における背景位相分布Ｐｅｄｇｅを算出する。
背景位相分布算出部３３５は、局所領域内の任意画素の位置を（ｘ、ｙ、ｚ）とし、式（６）を満たす多項式係数Ａ’ _{ｋ、ｌ、ｎ} を算出する。

Ｋはｘ方向の多項式における次数、Ｌはｙ方向の多項式における次数、Ｎはｚ方向の多項式における次数をそれぞれ表す。
局所領域内の中心画素の位置を（ｘ_０、ｙ_０、ｚ_０）とすると、対象組織の所定領域における背景位相分布Ｐｅｄｇｅは、式（７）により算出できる。

背景位相分布算出部３３５は、対象組織の所定領域における画素の全てに対して、式（７）を実施することにより、背景位相分布Ｐｅｄｇｅを算出する。

ステップＳ１１０６において、全領域局所位相分布算出部３３６は、総位相分布Ｐ_{ｔｏｔａｌ}と、対象組織の所定領域における背景位相分布Ｐｅｄｇｅとを組み合わせて、対象組織の全領域における局所位相分布Ｐ’_{ｌｏｃａｌ}を算出する。以下、対象組織の全領域における局所位相分布Ｐ’_{ｌｏｃａｌ}の算出法について説明する。

全領域局所位相分布算出部３３６は、制約付き最小二乗処理により、対象組織の全領域における局所位相分布の推定値Ｐ’_{ｌｏｃａｌ}を算出する。

Ｗ_１は脳表領域における重みが小さい重み画像、Ｗ_２は脳表領域における重みが大きい重み画像をそれぞれ表す。
式（８）における第一項は、ＳＨＡＲＰ法と同様、脳内部の領域における、平均値の性質に基づく制約項を表している。平均値の性質における条件を満たす制約項に対しては、脳内部の領域における重みを大きくし、脳表領域における重みを小さくする。
式（８）における第二項は、脳表領域における背景位相分布に基づく制約項を表している。脳表領域における背景位相分布に基づく制約項に対しては、脳表領域における重みを大きくし、脳内部の領域における重みを小さくする。

ここで、図７を参照して、重み画像Ｗ_１及び重み画像Ｗ_２と式（８）との関係について説明する。図７（ａ）は、重み画像Ｗ_１の一次元プロファイルを示している。図７（ｂ）は、重み画像Ｗ_２の一次元プロファイルを示している。横軸は位置を示し、縦軸は重みを示している。

図７において、領域７０１と領域７０２とを合わせた領域は、脳表領域を示し、領域７０３は、脳表領域を除いた脳内部の領域（脳表領域以外の領域）を示している。同様に、領域７０４と領域７０５とを合わせた領域は、脳表領域を示し、領域７０６は、脳表領域を除いた脳内部の領域（脳表領域以外の領域）を示している。なお、領域７０１（領域７０４）及び領域７０２（領域７０５）は、球状カーネルρ_ｒの半径ｒよりも大きくなるように設定される。

図７に示すように、全領域局所位相分布算出部３３６は、領域７０１での重みが０となるように設定し、領域７０３での重みが１となるように設定し、領域７０１から領域７０３へと遷移する領域７０２での重みが０から１の間の値（線形に変化する値）となるように設定する。また、全領域局所位相分布算出部３３６は、領域７０４での重みが１となるように設定し、領域７０６での重みが０となるように設定し、領域７０４から領域７０６へと遷移する領域７０５での重みが０から１の間の値（線形に変化する値）となるように設定する。

即ち、全領域局所位相分布算出部３３６は、重み画像Ｗ_１及び重み画像Ｗ_２を用いて、脳内部の領域における重みを大きくし、脳表領域における重みを小さくすることで、式（８）に示す平均値の性質に基づく制約項を、ＳＨＡＲＰ法と同等の精度で算出することができる（式（５）参照）。

また、全領域局所位相分布算出部３３６は、重み画像Ｗ_１及び重み画像Ｗ_２を用いて、脳表領域における重みを大きくし、脳内部の領域における重みを小さくすることで、式（８）に示す脳表領域における背景位相分布に基づく制約項を、局所多項式近似処理により算出された脳表領域における背景位相分布に近づけることができる。

更に、全領域局所位相分布算出部３３６は、重み画像Ｗ_１及び重み画像Ｗ_２を用いて、領域７０２及び領域７０６における重みを、０から１の間で線形に変化させることで、脳表領域（領域７０１と領域７０２とを合わせた領域、領域７０４と領域７０５とを合わせた領域）と、脳内部の領域（領域７０３、領域７０６）とを、滑らかに接続することができる。

なお、全領域局所位相分布算出部３３６は、式（８）に示すような制約付き最小二乗処理のみならず、式（９）に示すような正則化付最小二乗処理（局所位相分布の正則化項を最小化する条件）を用いて、対象組織の全領域（脳表領域を含めた脳全体）における局所位相分布を算出してもよい。

例えば、全領域局所位相分布算出部３３６は、局所位相分布Ｐ_{ｌｏｃａｌ}の画素値のばらつきを考慮した正則化項を挿入して、式（９）を実施することにより、対象組織の全領域における局所位相分布を算出することも可能である。

λ_１は、対象組織の所定領域における制約の大きさを調整するための正則化パラメータである。λ_２は、局所位相分布Ｐ_{ｌｏｃａｌ}のばらつきを、ある程度許容するための正則化パラメータである。
式（９）における第一項は、ＳＨＡＲＰ法と同様、脳内部の領域における、平均値の性質に基づく制約項を表している。
式（９）における第二項は、脳表領域における背景位相分布に基づく制約項を表している。
式（９）における第三項は、位相分布から局所位相変化成分を過剰に除去してしまうことを防ぐ役割を持つ正則化項を表している。

式（８）及び式（９）より、本実施形態に係る計算機１０９の演算処理によれば、脳内部の領域における局所位相分布の算出精度は、従来のＳＨＡＲＰ法と略同様の算出精度を維持しつつ、脳表領域における局所位相分布の算出精度を、従来と比較して格段に向上させることができることが示された。即ち、本実施形態に係る計算機１０９の演算処理によれば、脳内部の領域における平均値の性質に基づく制約項及び脳表領域における背景位相分布に基づく制約項から、脳表領域を含めた脳全体の局所位相分布を高精度に算出できる
ことが示された。

［所望分布算出：Ｓ１００４］
Ｓ１００４において、計算機１０９は所望分布算出処理を行う。

所望分布算出部３４０は、ステップＳ１００３で算出された局所位相分布Ｐ_{ｌｏｃａｌ}を用いて、ＱＳＭ法やＳＷＩ法などにより、定量的磁化率分布、磁化率強調画像などの所望分布を算出する。例えば、所望分布算出部３４０は、局所磁場算出部と磁化率分布算出部とを備え、ＱＳＭ法により、定量的磁化率分布を算出する。また、例えば、所望分布算出部３４０は、マスク作成部と磁化率強調画像作成部とを備え、ＳＷＩ法により、磁化率強調画像を算出する。以下、所望分布算出部３４０が行う具体的な処理例として、ＱＳＭ法を用いた定量的磁化率分布の算出法と、ＳＷＩ法を用いた磁化率強調画像の算出法について説明する。

≪ＱＳＭ法を用いた定量的磁化率分布の算出法≫
まず、ＱＳＭ法による定量的磁化率分布の算出法について説明する。ＱＳＭ法は、局所位相分布Ｐ_{ｌｏｃａｌ}が、生体組織間の磁化率差を反映することを利用して、局所位相分布Ｐ_{ｌｏｃａｌ}から、局所的な磁場分布を算出し、磁場と磁化率の関係式を用いて定量的磁化率分布を算出する手法である。

位置ベクトルをｒとすると、生体組織間の磁化率差によって生じる相対的な磁場分布は、式（１０）で表される。

gはプロトンの核磁気回転比、Ｂ０は静磁場強度をそれぞれ表す。
ここで、局所的な磁場分布は、静磁場に関するマクスウェル方程式より、式（１１）で表される。

ここで、d（ｒ）は局所磁場分布、c（ｒ）は生体内の磁化率分布、aはベクトル（ｒ’−ｒ）と静磁場方向との為す角度、ｄ（ｒ）は点ダイポール磁場をそれぞれ表す。

即ち、式（１１）に示すように、局所磁場分布d（ｒ）は、磁化率分布c（ｒ）と点ダイポール磁場ｄ（ｒ）との畳み込み積分で表すことができる。従って、式（１１）の両辺をフーリエ変換することにより、式（１１）は、式（１２）に変換できる。

ｋ＝（ｋ_ｘ、ｋ_ｙ、ｋ_ｚ）は、ｋ空間上の位置ベクトルを表す。
D（ｋ）は、局所磁場分布d（ｒ）のフーリエ成分を表す。
Ｄ（ｋ）は、点ダイポール磁場ｄ（ｒ）のフーリエ成分を表す。
Ｘ（ｋ）は、生体内の磁化率分布c（ｒ）のフーリエ成分を表す。

式（１２）に示すように、生体内の磁化率分布c（ｒ）のフーリエ成分Ｘ（ｋ）は、局所磁場分布d（ｒ）のフーリエ成分D（ｋ）を、点ダイポール磁場ｄ（ｒ）のフーリエ成分Ｄ（ｋ）で除算することによって算出できる。

しかしながら、式（１２）は、Ｄ（ｋ）＝０近傍の領域において、その逆数が発散してしまうため、直接的にＸ（ｋ）を算出することができない。このＤ（ｋ）＝０近傍の領域はマジックアングル領域と呼ばれ、静磁場方向に対しておよそ５４．７度の２倍の頂角を持つ逆双円錐領域となる。

ＱＳＭ法は、マジックアングル領域の存在により不良条件逆問題（ill-conditioned inverse problem）に帰着されるため、不良条件逆問題を解くための幾つかの解法が提案されている。不良条件逆問題を解く代表的な方法としては、例えば、磁場と磁化率との関係式に基づく制約条件下で、磁場分布から算出した磁化率分布に対して平滑化処理を行うことを繰り返す方法、磁場分布及び点ダイポール磁場のｋ空間上の演算により磁化率分布を算出するＴＫＤ（Truncated-based K-space Division）法、ＴＫＤ法で算出した磁化率分布と閾値処理により微細構造を抽出した磁化率分布とを繰り返し演算により合成するIterative ＳＷＩＭ（Susceptibility Weighted Imaging and Mapping）法、正則化付最小二乗法を用いたＭＥＤＩ（Morphology enabled dipole inversion）法、などが挙げられる。なお、所望分布算出部３４０は、これらの何れの方法を用いてもよい。

≪ＳＷＩ法を用いた磁化率強調画像の算出法≫
次に、ＳＷＩ法による磁化率強調画像の算出法について説明する。ＳＷＩ法は、局所位相分布Ｐ_{ｌｏｃａｌ}から、磁化率分布を強調する磁化率強調マスクを作成し、計測した強度画像（絶対値画像）に、磁化率強調マスクを乗算することで、磁化率強調画像を算出する手法である。

所望分布算出部３４０は、所定の磁化率強調マスクを作成し、この磁化率強調マスクを、所定回数乗算する。その後、所望分布算出部３４０は、任意のエコー時間の絶対値画像に、当該磁化率強調マスクを乗算することで、磁化率強調画像を算出する。磁化率強調画像とは、磁化率の高い所が暗くなり、磁化率の低い所が明るくなる画像である。操作者は、当該画像に基づいて、磁化率が高くなっているか、或いは低くなっているかを判断することができる。なお、所望分布算出部３４０は、公知のさまざまな手法を用いて、磁化率強調画像を算出することが可能であり、ＳＷＩ法に限定されるものではない。

［画像表示：Ｓ１００５］
Ｓ１００５において、計算機１０９は画像表示処理を行う。

計算機１０９は、ステップＳ１００４で算出された定量的磁化率分布や磁化率強調画像を、表示部１１０のディスプレイ画面上に表示させる。また、計算機１０９は、ステップＳ１００４で算出された各種データを外部記憶装置１１１或いは記憶装置に記憶させ、所望の表示装置に、定量的磁化率分布や磁化率強調画像を表示させることも可能である。

ここで、図８を参照して、ＭＲＩ装置１００の表示部１１０に表示される画像の一例について説明する。

図８（ａ）は、本実施形態に係る計算機１０９によって演算処理が施された画像（磁化率強調画像）であり、図８（ｂ）は、図８（ａ）の領域８０１における拡大図である。図８（ｃ）は、従来の計算機によって演算処理が施された画像（磁化率強調画像）であり、図８（ｄ）は、図８（ｃ）の領域８０２における拡大図である。

図８（ａ）及び図８（ｂ）に示す画像と図８（ｃ）及び図８（ｄ）に示す画像とを比較すると、図８（ａ）及び図８（ｂ）に示す画像では、脳表領域が鮮明に撮像されている一方で、図８（ｃ）及び図８（ｄ）に示す画像では、脳表領域が太い実線で囲まれたように撮像されていることがわかる。即ち、本実施形態に係る計算機１０９は、従来の計算機と比較して、太い実線で囲まれた領域の分だけ、より多くの情報を算出できることがわかる。

従って、操作者は、脳全体における生体組織の構造や組成などを診断する場合、図８（ｃ）及び図８（ｄ）に示す画像を用いるよりも、図８（ａ）及び図８（ｂ）に示す画像を用いることで、より適切な診断を行うことができる。更に、操作者は、くも膜下出血や脳梗塞など、脳の表面に発生した病気を診断する場合、図８（ａ）及び図８（ｂ）に示す画像を用いることで、正確な診断を行うことができる。即ち、本実施形態に係る計算機１０９が所定の演算処理を行うことで、操作者は、脳表領域において信頼性の高い画像を取得することができる。

以上、上述の各処理を行うことで、本実施形態に係る計算機１０９が行う演算処理が終了する。なお、上述の図４及び図６に示すフローは一例であり、その一部を省略することや、別の処理を追加することが可能であることは勿論である。

本実施形態によれば、計算機１０９が、複素画像から、例えば脳領域を抽出し、総位相分布を脳表領域と脳内部の領域とに分割して、脳表領域における背景位相分布を算出し、総位相分布と該背景位相分布とを組み合わせることで、脳全体の局所位相分布を算出する。これにより、脳表領域を含めた脳全体の局所位相分布を高精度に算出するＭＲＩ装置を実現することができる。

更に、本実施形態によれば、高精度な局所位相分布を用いて、定量的磁化率分布や磁化率強調画像を算出することができる。これにより、操作者は、鮮明な領域が拡張された画像を取得することができるため、ＭＲＩ装置の信頼性を高められる。

１００，１２０，１３０ 磁気共鳴イメージング装置
１０１ 被検体
１０２ 静磁場生成磁石
１０３ 傾斜磁場印加部
１０５ 送信部
１０６ 受信部
１０９ 計算機
１１０ 表示部
３３１ 対象組織抽出部
３３２ 位相分布算出部
３３３ 位相折り返し補正部
３３４ 領域分割部
３３５ 背景位相分布算出部
３３６ 全領域局所位相分布算出部（局所位相分布算出部）

Claims (6)

1. 被検体が配置される空間に静磁場を生成する静磁場生成磁石と、
   前記被検体に高周波磁場を送信する送信部と、
   前記高周波磁場の送信により前記被検体から発生する核磁気共鳴信号を受信する受信部と、
   前記核磁気共鳴信号に位置情報を付加する傾斜磁場を印加する傾斜磁場印加部と、
   前記送信部、前記受信部、及び前記傾斜磁場印加部を制御するとともに、前記核磁気共鳴信号に対して演算処理を行う計算機と、
   前記計算機により演算処理された画像を表示する表示部と、を備え、
   前記計算機は、
   前記核磁気共鳴信号によって構成された複素画像から前記被検体の対象組織を抽出した対象組織抽出画像を算出する対象組織抽出部と、
   前記複素画像から位相分布を算出する位相分布算出部と、
   前記位相分布の位相折り返しを補正する位相折り返し補正部と、
   前記対象組織抽出画像に基づいて、補正後の位相分布を、前記対象組織の所定領域と前記所定領域以外の領域とに分割する領域分割部と、
   前記所定領域の背景位相分布を算出する背景位相分布算出部と、
   前記補正後の位相分布と前記背景位相分布とを組み合わせて、前記対象組織の所定領域における局所位相分布を算出する局所位相分布算出部と、
   を備えることを特徴とする磁気共鳴イメージング装置。
2. 前記局所位相分布算出部は、
   任意画素の画素値と、前記任意画素を中心とした所定の半径を有する球内における画素の画素値の平均とが略等しいという条件に基づいて、前記局所位相分布を算出する、
   ことを特徴とする請求項１に記載の磁気共鳴イメージング装置。
3. 前記局所位相分布算出部は、
   前記局所位相分布の正則化項を最小化する条件に基づいて、前記局所位相分布を算出する、
   ことを特徴とする請求項２に記載の磁気共鳴イメージング装置。
4. 前記計算機は、
   前記局所位相分布を用いて、局所磁場分布を算出し、
   前記局所磁場分布及び磁場と磁化率との関係式を用いて、磁化率分布を算出する、
   ことを特徴とする請求項１から請求項３のいずれか一項に記載の磁気共鳴イメージング装置。
5. 前記計算機は、
   前記局所位相分布を用いて、前記磁化率分布を強調する磁化率強調マスクを作成し、
   前記複素画像の絶対値成分に、前記磁化率強調マスクを乗じて磁化率強調画像を算出する、
   ことを特徴とする請求項４に記載の磁気共鳴イメージング装置。
6. 前記領域分割部は、
   前記補正後の位相分布を、脳表領域と脳内部の領域とに分割する、
   ことを特徴とする請求項１から請求項５のいずれか一項に記載の磁気共鳴イメージング装置。

Images