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

Abstract

【課題】ＭＲＩを用いて、生体組織間の磁化率差に起因した局所磁場分布を算出するにあたり、演算時間が短く，ＳＮＲが高い局所周波数分布を算出する方法を提供する。【解決手段】ＭＲＩを用いて少なくとも２つ以上の異なるエコー時間で計測したマルチエコー複素画像を低解像度画像に変換する。低解像度マルチエコー複素画像の位相分布から、大域的な磁場変化に起因する大域的周波数分布と、受信や送信位相などを含むオフセット位相分布を分離する。算出した大域的周波数分布とオフセット位相分布を高解像度化する。計測したマルチエコー複素画像と高解像大域的周波数分布と高解像オフセット位相分布とから各エコーの局所周波数分布を算出する。各エコーの局所周波数分布を加重平均して最終的な局所周波数分布を算出する。【選択図】図６

Description

本発明は、磁気共鳴イメージング（ＭＲＩ：Ｍａｇｎｅｔｉｃ Ｒｅｓｏｎａｎｃｅ Ｉｍａｇｉｎｇ、以下、ＭＲＩという）技術に関する。特に、２つ以上の異なる複数のエコー時間で取得したマルチエコー複素画像を用いて生体組織の磁化率差に起因した局所磁場変化の分布を算出する画像処理技術に関する。

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

近年、ＭＲＩで測定できる物性値として、生体組織間の磁化率差が注目されている。磁化率とは、静磁場中の物質が磁気分極（磁化）する度合いを表した物性値である。生体内には、静脈血中のデオキシヘモグロビンや鉄たんぱく質などの常磁性体と、生体組織の大部分を占める水や石灰化の基となるカルシウムなどの反磁性体とがある。生体組織間の磁化率差を定量的に画像化したり、磁化率差を強調したりすることで、脳虚血疾患の診断やがんの放射線治療効果の予測、神経変性疾患の鑑別に適用できる可能性がある。

ＭＲＩを用いて生体組織間の磁化率差を画像化する方法は、定量的磁化率マッピング（ＱＳＭ：Ｑｕａｎｔｉｔａｔｉｖｅ Ｓｕｓｃｅｐｔｉｂｉｌｉｔｙ Ｍａｐｐｉｎｇ）法と呼ばれる。また、生体組織間の磁化率差を強調して画像化する方法は、磁化率強調画像化（ＳＷＩ：Ｓｕｓｃｅｐｔｉｂｉｌｉｔｙ Ｗｅｉｇｈｔｅｄ Ｉｍａｇｉｎｇ）法と呼ばれる。ＱＳＭ法は、計測したＭＲ画像の位相情報から生体組織間の磁化率差によって生じる局所的な磁場変化を算出し、磁場と磁化率の関係式を用いて定量的磁化率分布を算出する手法である。また、ＳＷＩは、局所的な磁場変化を強調する強調マスク画像を算出し、計測した強度画像（絶対値画像）にこの強度マスクを乗算することで磁化率が強調された画像を算出する手法である。

ＱＳＭ法やＳＷＩ法を用いて定量的磁化率分布や磁化率強調画像を算出するには、生体組織間の磁化率差によって生じる局所的な磁場変化を算出する必要がある。通常は、Ｇｒａｄｉｅｎｔ ｅｃｈｏ（ＧｒＥ）法を用いて、一つのエコー時間（ＴＥ）で計測した位相分布から、局所的な磁場変化の分布（局所磁場分布）を算出する。具体的には、計測した複素画像から、位相分布を算出し、算出した位相分布に生じている−π〜＋πの範囲に折り返された位相を除去する位相折り返し除去処理を実施する。その後、被検体の形状などに起因して生じる大域的磁場（背景磁場）変化を除去する背景除去処理を実施することで生体組織の磁化率差に起因した局所磁場を算出する。

局所磁場分布の画質は、計測した位相分布の信号対雑音比（ＳＮＲ：Ｓｉｇｎａｌ Ｎｏｉｓｅ Ｒａｔｉｏ）に依存する。ＭＲＩの位相分布は、対象とする生体組織の見かけの横磁化緩和時間Ｔ_２ ^＊に一致したＴＥで計測することで位相分布のＳＮＲが最大となることが知られている。一方、位相折り返し除去処理において、−π〜＋πで折り返された位相の折り返し除去を失敗なく実施するためにはＴＥを長くすることはできない。そのため一つのＴＥでＳＮＲが最適な位相分布を取得することはできない。そこで、複数のＴＥの画像を取得するマルチエコー計測を実施することで位相折り返し除去の失敗なく、ＳＮＲが最適な局所磁場分布を得られる。

マルチエコー計測によって得られた複数のＴＥの位相分布から局所磁場分布を算出する方法はいくつか提案されている。

代表的な方法として以下の二つの方法がある。

一つの方法（従来法１という）は、複数のＴＥで計測したマルチエコー画像の位相分布において、時間方向の位相折り返し除去処理を実施する。そして、時間方向の位相が線形変化することを利用して線形フィッティングすることで、静磁場不均一に起因した周波数成分を算出する。次いで、空間的な周波数折り返し除去処理を行い、その後、背景除去処理により局所磁場分布を算出する（例えば、非特許文献１参照）。

もう一つの方法（従来法２という）は、複数のＴＥで計測したマルチエコー画像の位相分布において、各エコーの位相分布に対して、個別に位相折り返し除去処理と背景除去処理を実施し、各エコーの局所磁場分布を算出する。そして、各エコーの局所磁場分布を加重平均することで最終的な局所磁場分布を算出する（例えば、特許文献１参照）。

Ｌｉｕ Ｔ他、 " Ｎｏｎｌｉｎｅａｒ Ｆｏｒｍｕｌａｔｉｏｎ ｏｆ ｔｈｅ Ｍａｇｎｅｔｉｃ Ｆｉｅｌｄ ｔｏ Ｓｏｕｒｃｅ Ｒｅｌａｔｉｏｎｓｈｉｐ ｆｏｒ Ｒｏｂｕｓｔ Ｑｕａｎｔｉｔａｔｉｖｅ Ｓｕｓｃｅｐｔｉｂｉｌｉｔｙ Ｍａｐｐｉｎｇ"、 Ｍａｇｎｅｔｉｃ Ｒｅｓｏｎａｎｃｅ ｉｎ Ｍｅｄｉｃｉｎｅ、２０１３年、６９巻、４６７−４７６頁

米国特許出願公開第２０１２／０３２１１６２号明細書

従来法１は、時間方向の位相折り返し除去処理を実施するため、位相の折り返しを精度よく除去できる利点がある。しかしながら、画素ごとに位相の線形フィッティング処理を実施する必要があるため、画素数に応じて演算時間が増大するという課題がある。また、ノイズの影響により、算出した局所磁場分布のＳＮＲが低いという課題がある。

一方、従来法２は、各エコーの局所磁場分布を加重平均するため、算出した局所磁場のＳＮＲが向上するという利点がある。しかしながら、各エコー画像で背景除去処理を実施する必要があるため、エコー数の増加に応じて演算時間が増大するという課題がある。また、長いＴＥで計測したエコー画像の場合、生体組織間の磁化率差が大きい領域では、位相折り返し除去が不正確となるという課題がある。

本発明は、上記事情に鑑みて成されたもので、ＭＲＩで計測した計測データを用いて、生体組織間の磁化率差に起因した局所磁場分布を算出するにあたり、演算時間が短く，ＳＮＲが高い局所周波数分布を算出する方法を提供することを目的とする。

本発明は、ＭＲＩを用いて少なくとも２つの異なるエコー時間で計測したマルチエコー複素画像を低解像度画像に変換する。低解像度マルチエコー複素画像の位相分布から、大域的な磁場変化に起因する大域的周波数分布と、受信や送信位相などを含むオフセット位相分布を分離する。算出した大域的周波数分布とオフセット位相分布を高解像度化する。計測したマルチエコー複素画像と高解像度大域的周波数分布と高解像度オフセット位相分布とから各エコーの局所周波数分布を算出する。各エコーの局所周波数分布を加重平均して最終的な局所周波数分布を算出する。

ＭＲＩを用いて、生体組織間の磁化率差に起因した局所磁場分布を算出するにあたり、短い演算時間で、ＳＮＲが高い局所周波数分布を得ることができる。

（ａ）〜（ｃ）は、それぞれ、本発明が適用されるＭＲＩ装置の外観図。 ＭＲＩ装置の構成の一実施形態を示す図 第一の実施形態の計算機の機能ブロック図 第一の実施形態の計算機による処理の流れを示す図 複素画像用データを取得するのに用いられるパルスシーケンスの一例を示す図 第一の実施形態における局所周波数分布算出処理の流れを示す図 第一の実施形態における分布画像算出部の機能ブロック図 第二の実施形態における位相分布分離処理の流れを示す図 第二の実施形態における位相分布分離処理部の機能ブロック図 第三の実施形態における位相分布分離処理の流れを示す図 第三の実施形態における位相分布分離処理部の機能ブロック図 第四の実施形態における位相分布分離処理の流れを示す図 第四の実施形態における位相分布分離処理部の機能ブロック図

本発明の実施形態の一例を説明する。以下、実施形態を説明するための全図において、特に断らない限り、同一機能を有するものは同一符号を付し、その繰り返しの説明は省略する。

＜第一実施形態＞
[ＭＲＩ装置の外観]
まず、本発明が適用されるＭＲＩ装置の実施形態について説明する。図１（ａ）〜図１（ｃ）は、ＭＲＩ装置の外観図である。図１（ａ）は、ソレノイドコイルで静磁場を生成するトンネル型磁石を用いた水平磁場方式のＭＲＩ装置１００である。図１（ｂ）は、開放感を高めるために磁石を上下に分離したハンバーガー型（オープン型）の垂直磁場方式のＭＲＩ装置１２０である。また、図１（ｃ）は、図１（ａ）と同じトンネル型磁石を用い、磁石の奥行を短くし、かつ、斜めに傾けることによって、開放感を高めたＭＲＩ装置１３０である。

本実施形態では、これらの外観を有するＭＲＩ装置のいずれを用いることもできる。但し、これらは一例であり、本実施形態のＭＲＩ装置はこれらの形態に限定されるものではない。本実施形態では、装置の形態やタイプを問わず、公知の各種のＭＲＩ装置を用いることができる。以下、特に区別する必要がない場合は、ＭＲＩ装置１００で代表する。

［ＭＲＩ装置の構成］
図２は、本実施形態のＭＲＩ装置１００の機能構成図である。本図に示すように、本実施形態のＭＲＩ装置１００は、被検体１０１が置かれる空間に静磁場を生成する、例えば、静磁場コイル、静磁場生成磁石などで構成される静磁場コイル１０２と、静磁場分布を調整するシムコイル１０４と、被検体１０１の計測領域に対し高周波磁場を送信する送信用高周波コイル１０５（以下、単に送信コイルという）と、被検体１０１から生じる核磁気共鳴信号を受信する受信用高周波コイル１０６（以下、単に受信コイルという）と、被検体１０１から生じる核磁気共鳴信号に位置情報を付加するために、ｘ方向、ｙ方向、ｚ方向それぞれに傾斜磁場を印加する傾斜磁場コイル１０３と、送信機１０７と、受信機１０８と、計算機１０９と、傾斜磁場用電源部１１２と、シム用電源部１１３と、シーケンス制御装置１１４と、を備える。

静磁場コイル１０２は、図１（ａ）、図１（ｂ）、図１（ｃ）にそれぞれ示した各ＭＲＩ装置１００、１２０、１３０の構造に応じて、種々の形態のものが採用される。

送信コイル１０５が照射する高周波磁場は、送信機１０７により生成される。すなわち送信コイル１０５と送信機１０７は送信部として機能する。また受信コイル１０６が検出した核磁気共鳴信号は、受信機１０８を通して計算機１０９に送られる。すなわち受信コイル１０６と受信機１８は受信部として機能する。なお、図２では、送信コイル１０５と受信コイル１０６とに別個のものを用いる場合を示しているが、送信コイル１０５と受信コイル１０６との機能を兼用する１つのコイルで構成してもよい。

傾斜磁場コイル１０３及びシムコイル１０４は、それぞれ傾斜磁場用電源部１１２及びシム用電源部１１３により駆動される。傾斜磁場コイル１０３と傾斜磁場用電源部１１２は傾斜磁場発生部として機能する。

シーケンス制御装置１１４は、傾斜磁場コイル１０３の駆動用電源である傾斜磁場用電源部１１２、シムコイル１０４の駆動用電源であるシム用電源部１１３、送信機１０７及び受信機１０８の動作を制御し、傾斜磁場、高周波磁場の印加および核磁気共鳴信号の受信のタイミングを制御する。制御のタイムチャートはパルスシーケンスと呼ばれ、計測に応じて予め設定され、後述する計算機１０９が備える記憶装置等に格納される。

計算機１０９は、ＭＲＩ装置１００全体の動作を制御するとともに、受信した核磁気共鳴信号に対して様々な演算処理を行う。本実施形態では、少なくとも２つ以上のエコー時間の複素画像や位相分布、磁場分布、磁化率分布、磁化率強調画像などを生成する。計算機１０９は、ＣＰＵ、メモリ、記憶装置などを備える情報処理装置であり、計算機１０９にはディスプレイ１１０、外部記憶装置１１１、入力装置１１５などが接続される。

ディスプレイ１１０は、演算処理で得られた結果等をオペレータに表示するインタフェースである。入力装置１１５は、本実施形態で実施する計測や演算処理に必要な条件、パラメータ等をオペレータが入力するためのインタフェースである。ユーザーは、入力装置１１５を介して、例えば、計測するエコーの数や、基準のエコー時間、エコー間隔などの計測パラメータを入力できる。外部記憶装置１１１は、計算機１０９内部の記憶装置とともに、計算機１０９が実行する各種の演算処理に用いられるデータ、演算処理により得られるデータ、入力された条件、パラメータ等を保持する。

［計算機の構成］
本実施形態のＭＲＩ装置では、計算機１０９が、複素画像を用いて、磁場分布、磁化率分布、磁化率強調画像などを生成するための演算を行う。このような演算を実現する計算機１０９の構成例を説明する。

図３は、本実施形態の計算機１０９の機能ブロック図である。本実施形態の計算機１０９は、高周波磁場パルスの照射に応じて被検体から発生する少なくとも２つ以上の異なるエコー時間の核磁気共鳴信号（エコー信号）を複素信号として計測する計測制御部３１０と、計測制御部３１０が計測した複素信号から画素値が複素数である複素画像を再構成する画像再構成部３２０と、画像再構成部３２０が再構成した複素画像から局所的周波数分布を算出する局所的周波数分布算出部３３０と、局所的周波数分布算出部３３０が算出した局所的周波数分布から、定量的磁化率分布や磁化率強調画像など目的とする分布画像を算出する分布画像算出部３４０とを備える。局所的周波数分布算出部３３０の詳細は後述する。

上述した計算機１０９の各部の機能は、記憶装置が保持するプログラム（ソフトウェア）を、ＣＰＵがメモリにロードして実行することにより実現される。各機能の処理に用いる各種のデータ、処理中に生成される各種のデータは、記憶装置あるいは外部記憶装置１１１に格納される。
また、計算機１０９が実現する各種の機能のうち、少なくとも一つの機能は、ＭＲＩ装置１００とは独立した、情報処理装置であって、ＭＲＩ装置１００とデータの送受信が可能な情報処理装置により実現されていてもよい。さらに、全部または一部の機能は、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（ｆｉｅｌｄ−ｐｒｏｇｒａｍｍａｂｌｅ ｇａｔｅ ａｒｒａｙ）などのハードウェアによって実現してもよい。

［計算機の処理］
上述した計算機１０９の構成を踏まえ、以下、計算機１０９が行う処理の概要を、図４を参照して説明する。

まず、計測制御部３１０は、予め定められたパルスシーケンスにしたがってシーケンス制御装置１１４を制御し、予め定めた少なくとも２つ以上の異なるエコー時間のエコー信号を計測する（ステップＳ１００１）。

その後、画像再構成部３２０は、得られたエコー信号を、ｋ空間上に配置してフーリエ変換することにより、複素画像Ｉ（ｎ）を再構成する（ステップＳ１００２）。これにより、異なるエコー時間毎の複素画像からなるマルチエコー複素画像が得られる。

局所的周波数分布算出部３３０は、複素画像から局所的周波数分布を算出する局所的周波数分布算出処理を行う（ステップＳ１００３）。

その後、分布画像算出部３４０は、算出した局所的周波数分布から、定量的磁化率分布や磁化率強調画像などの所望の分布画像の算出処理を行い（ステップＳ１００４）、算出した分布画像をディスプレイ１１０に表示する（ステップＳ１００５）。

なお、定量的磁化率分布や磁化率強調画像などの分布画像をディスプレイ１１０に表示する際、必要に応じて、ステップＳ１００３で算出した局所的周波数分布などのほか、所望分布算出処理の過程によって算出される画像をディスプレイ１１０に表示してもよい。

次に、各処理の詳細を説明する。
［計測：Ｓ１００１］
計測制御部３１０は、入力装置１１５を介してユーザーが入力したパラメータに基づいて設定されるパルスシーケンスに従って、シーケンス制御装置１１４を動作させ、予め定めたエコー時間（ＴＥ）の核磁気共鳴信号（エコー信号）を得る、計測を行う。シーケンス制御装置１１４は、計測制御部３１０からの指示に従って、ＭＲＩ装置１００の各部を制御して計測を行う。本実施形態では、少なくとも２つ以上のＮ個のエコー時間のエコー信号を得る。

ここで、計測制御部３１０が計測に用いるパルスシーケンスの例を説明する。本実施形態では、例えば、ＧｒＥ系のパルスシーケンスを用いる。このＧｒＥ系のパルスシーケンスは、生体組織間の磁化率差によって生じる局所的な磁場変化に鋭敏である。

図５に、ＧｒＥ系のパルスシーケンスの一例として、ＲＳＳＧ（ＲＦ−ｓｐｏｉｌｅｄ−Ｓｔｅａｄｙ−ｓｔａｔｅ Ａｃｑｕｉｓｉｔｉｏｎ ｗｉｔｈ Ｒｅｗｏｕｎｄ Ｇｒａｄｉｅｎｔ−Ｅｃｈｏ）−Ｍｕｌｔｉｅｃｈｏシーケンス５５０例を示す。本図において、ＲＦ、Ｇｓ、Ｇｐ、Ｇｒはそれぞれ、高周波磁場、スライス傾斜磁場、位相エンコード傾斜磁場、リードアウト傾斜磁場を表す。

ＲＳＳＧシーケンス５５０では、スライス傾斜磁場パルス５０１の印加とともに高周波磁場（ＲＦ）パルス５０２を照射し、被検体１０１内の所定のスライスの磁化を励起する。次いで磁化の位相にスライス方向および位相エンコード方向の位置情報を付加するためのスライスエンコード傾斜磁場パルス５０３および位相エンコード傾斜磁場パルス５０４を印加する。

画素内の核磁化の位相を分散させるディフェーズ用のリードアウト傾斜磁場パルス５０５を印加した後、リードアウト方向の位置情報を付加するためのリードアウト傾斜磁場パルス５０６、５０７、５０８、５０９を印加しながら核磁気共鳴信号（エコー）５１０，５１１，５１２，５１３をそれぞれ計測する。そして最後に、スライスエンコード傾斜磁場パルス５０３および位相エンコード傾斜磁場パルス５０４によってディフェーズされた核磁化の位相を収束させるリフェーズ用のスライスエンコード傾斜磁場パルス５１４および位相エンコード傾斜磁場パルス５１５を印加する。

計測制御部３１０は、以上の手順を、スライスエンコード傾斜磁場パルス５０３、５１４（スライスエンコード数ｋｓ）および位相エンコード傾斜磁場パルス５０４、５１５（位相エンコード数ｋｐ）の強度と、ＲＦパルス５０２の位相とを変化させながら、繰り返し時間ＴＲで繰り返し実行し、エコー時間毎に１枚の画像を得るために必要なエコーを計測する。なお、このとき、ＲＦパルス５０２の位相は、例えば、１１７度ずつ増加させる。また、図５において、ハイフン以下の数字は、繰り返しの何回目であるかを示す。

なお、計測される各エコーにおいて、血流などの流れの影響を補償するＦｌｏｗ Ｃｏｍｐｅｎｓａｔｉｏｎ傾斜磁場パルスを各軸に印加してもよい。

計測された各エコーはｋｒ、ｋｐ、ｋｓを座標軸とする３次元のｋ空間上（メモリ空間）に配置される。このとき、一つのエコーはｋ空間上でｋｒ軸に平行な１ラインを占める。このＲＳＳＧシーケンス５５０により得られる絶対値画像は、ＴＥが短いエコーの場合はＴ１（縦緩和時間）強調画像、ＴＥが長いエコーの場合は画素内の位相分散を反映したＴ２＊強調画像となる。

ここでは、ｋ空間の座標軸に平行にデータを取得するカーテシアン撮像の一つであるＲＳＳＧシーケンス５５０を例示するが、本実施形態では、任意のシーケンスを用い、任意のｋ空間領域のデータを取得することができる。例えば、ｋ空間において回転状にデータを取得するラジアルスキャンなど、ノンカーテシアン撮像を用いてもよい。また複数のＴＥのエコープラナー型のｋ空間走査法（マルチエコーエコープラナーイメージング法）でもよい。

さらに、複数のＴＥで取得した別々の複素画像や位相分布から一枚の複素画像や位相分布を作成するシーケンスを用いてもよい。例えば、複数のＴＥで得た信号から画像を得るＤｉｘｏｎ法を実現可能なパルスシーケンスを用いる場合、水と脂肪とを分離する過程で静磁場不均一周波数分布およびオフセット位相分布が得られる。

［画像再構成：Ｓ１００２］
次に、画像再構成部３２０は、ステップＳ１００１で計測した少なくとも２つの異なるエコー時間ＴＥ_ｎ（ｎはエコー番号を表し、１以上の整数）のエコー信号をｋ空間上に配置し、フーリエ変換する。これにより、画像再構成部３２０は、ＴＥ毎に、各画素値が複素数である複素画像Ｉ（ｎ）を算出する。

［局所的周波数分布算出処理：Ｓ１００３］
次に、局所的周波数分布算出部３３０は、画像再構成部３２０が再構成した複素画像Ｉ（ｎ）から、局所的周波数分布を算出する。局所的周波数分布は、生体組織間の磁化率差によって生じる局所的な周波数変化を捉えたものであり、ＱＳＭ法やＳＷＩ法を適用し、定量的磁化率分布や、磁化率強調画像を算出するために必要となる。

局所的周波数分布の具体的な算出方法を説明する前に、まず複素画像と局所的周波数分布との関係について説明する。

上述の通り複素画像は、異なるＴＥで計測したエコーから、ＴＥ毎に画像を再構成したものであり、絶対値成分と位相成分（位相分布）を含む。本実施形態の局所的周波数分布算出処理では、位相分布を用いて局所的周波数分布を算出する。

ここで、ＴＥの数をＮとすると、ｎ番目のＴＥで計測した複素画像複素画像Ｉ（ｎ）の位相分布Ｐ（ｎ）は、以下の式（１）に示す各成分に分解できる。

式中、Ｐ_{ｉｎｈｏｍｏ}は、生体組織間の磁化率差や被検体形状などによって生じる静磁場不均一に起因する位相分布、Ｐ_{ｏｆｆｓｅｔ}は計測時に生じる受信位相や送信位相などのオフセット位相分布をそれぞれ表す。

静磁場不均一に起因する位相分布Ｐ_{ｉｎｈｏｍｏ}には、生体組織間の磁化率差によって生じる局所磁場変化に起因した局所的位相分布と、被検体形状などよって生じる大域的な磁場変化に起因にした大域的位相分布（背景位相分布とも呼ばれる）が含まれる。従ってＰ_{ｉｎｈｏｍｏ}は式（２）に示す各成分に分解できる。

式中、Ｐ_{ｌｏｃａｌ}は、生体組織間の磁化率差によって生じる局所磁場変化に起因した局所的位相分布、Ｐ_{ｇｌｏｂａｌ}は、被検体形状などよって生じる大域的な磁場変化に起因にした大域的位相分布をそれぞれ表す。

一方、複素画像における位相Ｐと周波数Ｆとには所定の関係（Ｐ＝２πＦ×ＴＥ）があり、局所的位相分布Ｐ_{ｌｏｃａｌ}と、大域的位相分布Ｐ_{ｇｌｏｂａｌ}は、ＴＥ_ｎ（ｎ番目のＴＥを示す）を用いて式（３）、（４）でそれぞれ表される。

式中、Ｆ_{ｌｏｃａｌ}は、生体組織間の磁化率差によって生じる局所磁場変化に起因した局所的周波数分布、Ｆ_{ｇｌｏｂａｌ}は、被検体形状などよって生じる大域的な磁場変化に起因にした大域的周波数分布をそれぞれ表す。なお、以下では、Ｆ_{ｌｏｃａｌ}とＦ_{ｇｌｏｂａｌ}の和をＦ_{ｉｎｈｏｍｏ}と表し、静磁場不均一周波数分布と呼ぶ。以上の式（１）〜（４）から、オフセット位相分布と大域的周波数分布を算出することにより、局所的周波数分布が求められることがわかる。

局所的周波数分布算出部３３０は、上述した式（１）〜（４）に基づき、ＧｒＥ法で計測したマルチエコー複素画像の位相分布Ｐ（ｎ）から、局所的周波数分布Ｆ_{ｌｏｃａｌ}を算出する。その際、ＳＮＲを低下させることなく局所的周波数分布算出の負担の低減を実現する。

本実施形態における局所的周波数分布算出部３３０の処理（Ｓ１００３）の詳細を、図６に示す。図示するように、局所的周波数分布算出部３３０は、まず、複素画像Ｉ（ｎ）をもとの解像度より低い解像度（低解像度）の複素画像ｉ（ｎ）に変換する（Ｓ１１０１）。そして、低解像複素画像ｉ（ｎ）から低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を分離する（Ｓ１１０２）。次いで、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を、計測した複素画像Ｉ（ｎ）と同じ解像度となるように高解像度化して大域的周波数分布Ｆ_{ｇｌｏｂａｌ}とオフセット位相分布Ｐ_{ｏｆｆｓｅｔ}を算出する（Ｓ１１０３）。その後、計測した複素画像Ｉ（ｎ）から大域的周波数分布Ｆ_{ｇｌｏｂａｌ}とオフセット位相分布Ｐ_{ｏｆｆｓｅｔ}を除去することで、各エコーの局所的周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）を算出する（Ｓ１１０４）。そして、各エコーの局所的周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）を加重平均することで最終的な一つの局所的周波数分布Ｆ_{ｌｏｃａｌ}を得る（Ｓ１１０５）。

これらの処理は、図３に示す局所周波数分布算出部３３０の各部により実現することができる。すなわち、図示するように、局所周波数分布算出部３３０は、複素画像Ｉ（ｎ）を低解像複素画像ｉ（ｎ）に変換する低解像度変換部（第一の解像度変換部）３３１と、低解像複素画像ｉ（ｎ）から低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を分離する位相分布分離部３３２と、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を、計測した複素画像Ｉ（ｎ）と同じ解像度となるように高解像度化して大域的周波数分布Ｆ_{ｇｌｏｂａｌ}とオフセット位相分布Ｐ_{ｏｆｆｓｅｔ}を算出する高解像度変換部（第二の解像度変換部）３３３と、計測した複素画像Ｉ（ｎ）から大域的周波数分布Ｆ_{ｇｌｏｂａｌ}とオフセット位相分布Ｐ_{ｏｆｆｓｅｔ}を除去することで、各エコーの局所的周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）を算出する位相除去部３３４と、各エコーの局所的周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）を加重平均することで最終的な一つの局所的周波数分布Ｆ_{ｌｏｃａｌ}を得る加重平均部３３５と、を備える。

以下、本実施形態の局所周波数分布算出部３３０の各部が行う処理内容を、図６を参照して説明する。

［低解像度化：Ｓ１１０１］
低解像度変換部３３１は、計測したマルチエコー複素画像Ｉ（ｎ）を線形補間やＣｕｂｉｃ補間など公知の手法を用いて、計測した複素画像Ｉ（ｎ）の分解能以下となる所望の解像度に変換し、低解像複素画像ｉ（ｎ）を算出する。なお、低解像変換部３３１は、計測した被検体の構造に併せて、被検体内部の分解能を粗く、その辺縁を細かくするなど、領域によって任意の解像度に変換した低解像複素画像ｉ（ｎ）を算出してもよい。あるいは、複数の解像度の複素画像ｉ（ｎ）を算出してもよい。

［位相分布分離：Ｓ１１０２］
位相分布分離部３３２は、低解像複素画像ｉ（ｎ）から低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を分離する。この分離処理の方法には、位相差を求めることでオフセット位相分布を除去した位相差画像を得て、この位相差画像から大域的周波数分布ｆ_{ｇｌｏｂａｌ}を分離し、さらに低解像複素画像ｉ（ｎ）と大域的周波数分布ｆ_{ｇｌｏｂａｌ}を用いてオフセットｐ_{ｏｆｆｓｅｔ}を分離する方法や、信号論理式を用いたフィッティングによって、一度に低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を算出する方法などを採りえる。これら分離処理の具体的手法は、他の実施形態で説明する。

［高解像度化：Ｓ１１０３］
高解像度変換部３３３は、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を、線形補間やＣｕｂｉｃ補間など公知の手法を用いて、計測した複素画像Ｉ（ｎ）と同じ解像度となるように高解像度化して大域的周波数分布Ｆ_{ｇｌｏｂａｌ}とオフセット位相分布Ｐ_{ｏｆｆｓｅｔ}を算出する。

なお、図６に点線で囲って示すように（高解像度化変換前処理Ｓ１１０３−１）、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を高解像度化する前に、各低解像度画像の辺縁を外挿してから高解像度化してもよい。すなわち低解像度画像では画像の辺縁の情報が欠落している場合があるので、そのまま高解像度にした場合、高解像度画像に不連続（断絶）を生じる可能性がある。辺縁の情報をそれに近接する画素の情報で外挿入することにより、このような不連続を防止できる。外挿の方法は、特に限定されないが、線形補間やＣｕｂｉｃ補間の他、多項式フィッティングなどの各種方法を用いることができる。

さらに、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を高解像度化する前に平滑化処理によりノイズを除去してから高解像度化してもよい（平滑部の機能）。平滑化の方法は、任意の窓サイズの平均値フィルタやガウシアンフィルタなど各種方法を用いてもよい。ｆ_{ｇｌｏｂａｌ}やｐ_{ｏｆｆｓｅｔ}は、通常、緩やかに変化するので、平滑化処理を行うことでＳＮＲを向上することができる。

［位相除去：Ｓ１１０４］
位相分布除去部３３４は、計測した複素画像Ｉ（ｎ）から大域的周波数分布Ｆ_{ｇｌｏｂａｌ}とオフセット位相分布Ｐ_{ｏｆｆｓｅｔ}を除去することで、各エコーの局所的周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）を算出する。各エコーの局所的周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）は、次の式（５）〜（７）に基いて算出できる。

まず計測した複素画像Ｉ（ｎ）は、各エコーの局所的周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）と、大域的周波数分布Ｆ_{ｇｌｏｂａｌ}と、オフセット位相分布Ｐ_{ｏｆｆｓｅｔ}を用いて式（５）で表される。

式中、ｊは虚数、Ｍ_０はプロトン密度、ｋは計測シーケンスに依存した係数、である。
したがって、計測した複素画像Ｉ（ｎ）から、大域的周波数分布Ｆ_{ｇｌｏｂａｌ}と、オフセット位相分布Ｐ_{ｏｆｆｓｅｔ}による位相変化を複素除算した複素画像Ｉ’（ｎ）は、式（６）で表される。

ここで、式（７）に示すように、複素画像Ｉ’（ｎ）の位相成分をとり、位相成分から周波数成分に変換することで各エコーの局所的周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）を算出することができる。

［加重平均：Ｓ１１０５］
加重平均部３３５は、Ｓ１１０４で算出した各エコーの局所的周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）を、式（８）に示すように加重平均することで最終的な一つの局所的周波数分布Ｆ_{ｌｏｃａｌ}を得る。

式中、ｗ（ｎ）は加重平均に用いる重みであり、加重平均部３３５、或いは、計算機１０９内の記憶装置に予め設定しておくことができる。重みｗ（ｎ）は、例えば、複素画像Ｉ（ｎ）の絶対値成分｜Ｉ（ｎ）｜とエコー時間ＴＥ_ｎを用いて式（９）で表すことができる。

なお、重みｗ（ｎ）は式（９）で表すものに限られない。絶対値成分｜Ｉ（ｎ）｜の代わりに、式（１０）に示すように、見かけの横磁化緩和速度Ｒ_２ ^＊（横緩和時間Ｔ_２ ^＊の逆数）を用いてもよい。

以上の手順により、局所的周波数分布Ｆ_{ｌｏｃａｌ}を算出することができる。即ち、図４のステップＳ１００３の処理（局所的周波数分布算出処理）が完了する。

以上の手順により、信号フィッティングを用いることなく大域的周波数分布ｆ_{ｇｌｏｂａｌ}とオフセット位相分布ｐ_{ｏｆｆｓｅｔ}を算出することができる。したがって、非特許文献１に記載されるような信号フィッティング（従来法１）に比べて演算時間を短縮できる。そして、複素画像毎に算出したものを加重平均することで、精度よく大域的周波数分布ｆ_{ｇｌｏｂａｌ}とオフセット位相分布ｐ_{ｏｆｆｓｅｔ}を算出できる。また、本実施形態によれば、演算時間のかかる大域的周波数分布算出処理が一回で済むため、特許文献１に記載されるような個別処理法（従来法２）に比べて演算時間を短縮できる。また、オフセット位相分布ｐ_{ｏｆｆｓｅｔ}分離した上で大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出しているため、大域的周波数分布ｆ_{ｇｌｏｂａｌ}の算出精度が向上する。

［分布画像算出：Ｓ１００４］
図４に戻り、分布画像算出部３４０が行う処理（Ｓ１００４）について説明する。分布画像は、例えば、磁場分布や、定量的磁化率分布、磁化率強調画像などを含み、Ｓ１００３で算出した局所的周波数分布Ｆ_{ｌｏｃａｌ}から、ＱＳＭ法やＳＷＩ法を用いて算出することができる。これら処理を行う分布画像算出部３４０の機能ブロック図を、図７に示す。

図示する例では、分布画像算出部３４０は、定量的磁化率分布算出部７１０と磁化率強調画像算出部７２０とを有しているが、これらの一方のみを有していてもよい。定量的磁化率分布算出部７１０は、局所的周波数分布算出部３３０が算出した局所的周波数分布を用いて、局所磁場分布を算出する局所磁場算出部７１１、及び、局所磁場算出部７１０が算出した局所磁場と、磁場と磁化率との関係式とを用いて磁化率分布を算出する磁化率分布算出部７１２を有する。磁化率強調画像算出部７２０は、局所的周波数分布算出部３３０が算出した局所的周波数分布を用いて、磁化率強調マスクを作成するマスク作成部７２１、及び、複数の複素画像の絶対値成分に、マスク生成部７２１が作成した磁化率強調マスクを乗じて磁化率強調画像を作成する磁化率強調画像作成部７２２を有する。

以下、分布画像算出部３４０が行う具体的な処理例として、ＱＳＭ法を用いた定量的磁化率分布の算出（定量的磁化率分布算出部７１０の処理）と、ＳＷＩ法を用いた磁化率強調画像の算出（磁化率強調画像算出部７２０の処理）とを説明する。

［ＱＳＭ法を用いた定量的磁化率分布の算出法］
ＱＳＭ法では、生体組織間の磁化率差によって生じる局所的な磁場変化を計算し、磁場と磁化率の関係式を用いて定量的磁化率分布を算出する。このため、まず局所磁場算出部７１１は、局所的周波数分布Ｆ_{ｌｏｃａｌ}を用いて磁場変化を算出する。

ここで上述のように算出した局所的周波数分布Ｆ_{ｌｏｃａｌ}を用いたとき、位置ベクトルをｒとすると組織間の磁化率差によって生じる相対的な磁場変化（磁場分布）δ（ｒ）は、以下の式（１１）で表される。

式中、γはプロトンの核磁気回転比、Ｂ_０は静磁場強度をそれぞれ表す。

次いで磁化率分布算出部７１２が、磁場分布δ（ｒ）（局所磁場）と、磁場と磁化率との関係式とを用いて磁化率分布を算出する。
ここで磁場分布δ（ｒ）は、静磁場に関するマクスウェル方程式より、生体内の磁化率分布χ（ｒ）を用いて以下の式（１２）で表される。

式中、αは、ベクトル（ｒ’−ｒ）と静磁場方向とのなす角度、ｄ（ｒ）は、点ダイポール磁場をそれぞれ表す。

式（１２）に示すように、磁場分布δ（ｒ）は、磁化率分布χ（ｒ）と点ダイポール磁場ｄ（ｒ）との畳み込み積分で表される。したがって、式（１２）の両辺をフーリエ変換することにより、式（１２）は以下の式（１３）に変換される。

ここで、ｋ＝（ｋ_ｘ、ｋ_ｙ、ｋ_ｚ）は、ｋ空間上の位置ベクトル、Δ（ｋ）、Ｘ（ｋ）、Ｄ（ｋ）、は、磁場分布δ（ｒ）、磁化率分布χ（ｒ）、点ダイポール磁場ｄ（ｒ）のフーリエ成分をそれぞれ表す。

式（１３）に示すように、磁化率分布のフーリエ成分Ｘ（ｋ）は、磁場分布のフーリエ成分Δ（ｋ）を点ダイポール磁場のフーリエ成分Ｄ（ｋ）で除算することによって算出できる。しかしながら、式（１３）は、Ｄ（ｋ）＝０近傍の領域において、その逆数が発散してしまうため、直接的にＸ（ｋ）を算出することができない。

このＤ（ｋ）＝０となる領域はマジックアングルと呼ばれ、磁場方向に対しておよそ５４．７°の２倍の頂角を持つ逆双円錐領域となる。マジックアングルの存在により磁場分布から磁化率分布を推定するＱＳＭ法は、不良条件逆問題（ｉｌｌ−ｃｏｎｄｉｔｉｏｎｅｄ ｉｎｖｅｒｓｅ ｐｒｏｂｌｅｍ）に帰着され、いくつかの解法が提案されている。

その代表的な方法として、磁場と磁化率との関係式に基づく制約条件下で、磁場分布から算出した磁化率分布に対して平滑化処理を行うことを繰り返す方法（本発明者らによる特願２０１４−２２８８４３号記載の方法）や、磁場分布および点ダイポール磁場のｋ空間上の演算により磁化率分布を算出するＴＫＤ（Ｔｒｕｎｃａｔｅｄ−ｂａｓｅｄ Ｋ−ｓｐａｃｅ Ｄｉｖｉｓｉｏｎ）法、ＴＫＤ法で算出した磁化率分布と、閾値処理により微細構造を抽出した磁化率分布とを繰り返し演算により合成するＩｔｅｒａｔｉｖｅ ＳＷＩＭ（Ｓｕｓｃｅｐｔｉｂｉｌｉｔｙ Ｗｅｉｇｈｔｅｄ Ｉｍａｇｉｎｇ ａｎｄ Ｍａｐｐｉｎｇ）法、正則化付き最小二乗法を用いたＭＥＤＩ法（Ｍｏｒｐｈｏｌｏｇｙ ｅｎａｂｌｅｄ ｄｉｐｏｌｅ ｉｎｖｅｒｓｉｏｎ）法がある。
本実施形態の磁化率強調画像算出部７２０は、これらの方法を用いて定量的磁化率分布を算出する。本実施形態では、いずれの方法を用いても定量的磁化率分布を算出してもよい。

［ＳＷＩ法を用いた磁化率強調画像の算出法］
以下に、磁化率強調画像算出部７２０で行われる、ＳＷＩ法による磁化率強調画像の算出法について説明する。
本実施形態では、まず、マスク生成部７２１が、局所的周波数分布Ｆ_{ｇｌｏｂａｌ}から磁化率を強調する磁化率強調マスクを作成する。具体的には、任意の周波数ｆ_ｍ以下の値を０、０以上の周波数を１、その間を線形に結んだ値に設定する強調マスクを作成する。

その後、磁化率強調画像作成部７２で、強調マスクを任意の回数乗算したのち、任意のＴＥの複素画像の絶対値成分（絶対値画像）に乗算することで、磁化率強調画像を算出する。
なお、本実施形態の磁化率強調画像の算出法は上記の方法に限られない。公知のさまざまな手法を適用してもよい。

［画像表示：Ｓ１００５］
分布画像算出部３４０が算出した定量的磁化率分布や磁化率強調画像は、ディスプレイ１１０（図１）に表示することができる。或いは外部記憶装置１１１に画像データとして格納し、所望の表示装置で表示されてもよい。

なお図４及び図６に示すフローは一例であり、その一部を省くことや、別な処理を追加することも本実施形態に含まれる。例えば、図４の画像表示Ｓ１００５を省くこともあり得るし、例えば、ＴＥの異なる複素画像が２つの場合など、図６の加重平均Ｓ１１０５を省くことも可能である。

以上、本実施形態のＭＲＩ装置の構成と、その計算機が行う諸処理について説明した。本実施形態のＭＲＩ装置は、静磁場内に配置された被検体に高周波磁場パルスを送信する送信部と、前記被検体が発生する核磁気共鳴信号を受信する受信部と、前記受信部が受信した核磁気共鳴信号を処理し、計測データを作成する信号処理部と、前記信号処理部が作成した計測データに演算を施す計算機と、を備え、計算機が以下の処理を実現する機能を備えることが特徴である。

すなわち、本実施形態の計算機は、少なくとも２つ以上の異なるエコー時間で計測したマルチエコー複素画像を低解像度画像に変換する。低解像度マルチエコー複素画像の位相分布から、大域的な磁場変化に起因する大域的周波数分布と、受信や送信位相などを含むオフセット位相分布を分離する。算出した大域的周波数分布とオフセット位相分布を高解像度化する。計測したマルチエコー複素画像と高解像大域的周波数分布と高解像オフセット位相分布とから各エコーの局所周波数分布を算出する。各エコーの局所周波数分布を加重平均して最終的な局所周波数分布を算出する。

このため、計算機は、ＭＲＩ装置で取得した複数の複素画像データを、それぞれ、当該画像データより低解像度に変換する第一の解像度変換部と、前記第一の解像度変換部で処理された低解像度画像から、大域的周波数分布とオフセット位相分布とを分離する位相分布分離部と、前記位相分布分離部によって分離された前記大域的周波数分布及びオフセット位相分布を前記複数の画像データと同じ解像度に変換する第二の解像度変換部と、前記第二の解像度変化部で処理された前記大域的周波数分布及びオフセット位相分布と、前記複数の画像データとを用いて、局所的周波数分布を算出する局所的周波数分布算出部とを備える。

また必須ではないが、本実施形態の計算機は、低解像度オフセット位相分布及び／又は低解像大域的周波数分布を高解像度に変換する際に、平滑化する手段や辺縁のデータを外挿する手段を備えることができる。

上述した機能を持つ計算機は、ＭＲＩ装置の一部であってもよいし、ＭＲＩ装置から独立した演算処理装置であってもよい。即ち本実施形態は、ＭＲＩ装置のみならず、上述した機能を持つ計算機（演算処理装置）も包含する。

本実施形態のＭＲＩ装置（或いは演算処理装置）及び画像処理方法によれば、簡易に高い精度で局所的周波数分布を算出することができ、それを用いて精度のよい定量的磁化率分布や磁化率強調画像などの所望分布画像を得ることができる。

以上説明した第一の実施形態を踏まえ、計算機の局所的周波数分布算出部、特に位相分布分離部の機能を具体化した実施形態について、説明する。なお第一実施形態の説明で参照した図１〜図６に表される構成や処理は、特に断らない限り、以下説明する各実施形態に共通しており、適宜、これらの図面を参照するが、これら図面に記載された各要素についての重複する説明は省略する。

＜第二の実施形態＞
本実施形態は、位相分布分離部３３２（図３）が、低解像複素画像ｉ（ｎ）から低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を分離する際、位相差画像を作成し、オフセット位相分布を排除した静磁場不均一周波数分布を求め、この静磁場不均一周波数分布から、大域的周波数分布ｆ_{ｇｌｏｂａｌ}を分離する。

その前提として、本実施形態においても、図６のフローチャートに示すように、低解像度化処理（Ｓ１１０１）が行われ、また位相分布分離処理（Ｓ１１０２）に続いて、高解像度化処理（Ｓ１１０３）、位相分布除去処理（Ｓ１１０４）等が行われる。図８に、本実施形態の位相分布分離処理（図６：Ｓ１１０２）の概要を示す。

まず、位相分布分離部３３２は、低解像度変換部３３１にて複素画像Ｉ（ｎ）を低解像度化した低解像複素画像ｉ（ｎ）において、隣り合うエコー間の位相差分布Δｐ（ｎ）を（Ｎ−１）セット算出する（Ｓ１２０１）。そして、算出した位相差分布Δｐ（ｎ）を加重平均して一つの位相差分布Δｐを算出する（Ｓ１２０２）。位相差分布Δｐに対して位相アンラップ処理を実施し（Ｓ１２０３）、周波数成分に換算し静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}を算出する（Ｓ１２０４）。周波数分布分離による大域的周波数分布算出処理により、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}から大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する（Ｓ１２０５）。低解像複素画像ｉ（ＴＥ_ｎ）から、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}による位相変化を除去することで、各エコーのオフセット位相分布ｐ_{ｏｆｆｓｅｔ}（ｎ）を算出する（Ｓ１２０６）。各エコーのオフセット位相分布ｐ_{ｏｆｆｓｅｔ}（ｎ）を加重平均することでオフセット位相分布ｐ_{ｏｆｆｓｅｔ}を算出する（Ｓ１２０７）。以上の手順により、信号フィッティングを用いずに低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像度オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を算出することができる。

これらを実現するための位相分布分離部３３２の機能ブロック図を図９に示す。本図に示すように、位相分布分離部３３２は、隣り合うエコー間の位相差分布Δｐ（ｎ）を（Ｎ−１）セット算出する位相差算出部９０１と、算出した位相差分布Δｐ（ｎ）を加重平均した位相差分布Δｐに対して位相アンラップ処理を実施する位相アンラップ部９０２と、位相差成分を周波数成分に換算し、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}を算出する周波数変換部９０３と、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}から、大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する大域的周波数分布算出部（周波数分布分離部）９０４と、低解像複素画像ｉ（ｎ）から、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}による位相変化を除去することで、各エコーのオフセット位相分布ｐ_{ｏｆｆｓｅｔ}（ｎ）を算出するオフセット位相算出部９０５と、を備える。

以下、位相分布分離部３３２の各部が行う処理内容を、図８を参照して詳述する。

［位相差分布算出：Ｓ１２０１］
位相差算出部９０１において、低解像度マルチエコー複素画像ｉ（ｎ）における隣り合うエコー間の位相差分布Δｐ（ｎ）を（Ｎ−１）セット算出する。ここで、（ｎ＋１）番目のエコー時間ＴＥ_ｎ＋１における低解像複素画像ｉ（ｎ＋１）と、ｎ番目のエコー時間ＴＥ_ｎにおける低解像複素画像ｉ（ｎ）は、低解像度の静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}と、低解像オフセット位相ｐ_{ｏｆｆｓｅｔ}を用いて、それぞれ式（１４）と式（１５）で表される。

ここで、ＴＥ_ｎ＋１とＴＥ_ｎとの間のエコー間隔（エコー時間差分）をΔＴＥ_ｎとするとき、ｉ（ｎ＋１）をｉ（ｎ）で複素除算した複素比画像Δｉ（ｎ）は式（１６）で表される。

したがって、隣り合うエコー間の位相差分布Δｐ（ｎ）は式（１７）より算出することができる。

式（１６）、（１７）にそれぞれ示すように、隣り合うエコー間の位相差分布Δｐ（ｎ）は、低解像度のオフセット位相分布ｐ_{ｏｆｆｓｅｔ}が除去され、静磁場不均一に起因した静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}による位相変化のみに依存することが分かる。
Ｓ１２０１では、この処理を（Ｎ−１）回繰り返すことで隣り合うエコー間の位相差分布Δｐ（ｎ）を（Ｎ−１）セット算出する。

［加重平均：Ｓ１２０２］
次に、加重平均部３３５において、式（１８）に示すように、隣り合うエコー間の位相差分布Δｐ（ｎ）を加重平均することで一つの位相差分布Δｐを算出する。

式中、ｗ（ｎ）は重みであり、第一実施形態のステップＳ１００３（図４）において、各エコーの局所周波数分布Ｆ_{ｌｏｃａｌ}（ｎ）を加重平均する際に用いた重み（式（９）や式（１０）に示す重み）など、任意の重みを用いてもよい。なお、加重平均部３３５における加重平均処理は、式（１９）に示すように複素画像上で加重平均したのちに、位相成分を算出する方法でもよい。

この加重平均処理により算出した位相差分布ΔｐのＳＮＲが改善するため、後述の処理の精度が向上する。

［位相アンラップ：Ｓ１２０３］
次に、位相アンラップ部９０２において、位相差分布Δｐに対し、−πからπの範囲を超えて折り返された位相を除去するための位相アンラップ処理を実行する。位相差分布Δｐにおける一部の領域では、−πからπの範囲を超えた位相値が−πからπの範囲内に折り返される。そのため、撮像部位（例えば、頭部など）の全領域において正確な位相を求めるためには、これらの折り返しを補正する必要がある。本実施形態では、例えば領域拡大法などを用い、−πからπの範囲に折り返されている位相値を補正する。

［周波数変換：Ｓ１２０４］
次に、周波数変換部９０３において、式（２０）に示すように、位相差分布Δｐを、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}に変換する。

［大域的周波数分布算出（周波数分布分離）：Ｓ１２０５］
次に、大域的周波数分布算出部９０４において、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}から、生体組織間の磁化率差などに起因して生じる局所的周波数変化と、生体形状などに起因して生じる大域的周波数変化とを分離し、大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する。

大域的周波数分布を算出する典型的な方法として、Ｓｐｈｅｒｉｃａｌ Ｍｅａｎ Ｖａｌｕｅ（ＳＭＶ）フィルタ処理により算出する方法があり、これを採用することができる。ＳＭＶフィルタ処理は、大域的周波数分布ｆ_{ｇｌｏｂａｌ}が球面調和関数展開できることを利用する。球面調和関数は、任意の点ｒにおける値が、その周囲の任意の半径を持つ球内の平均値に等しい、という特性を持つ。ＳＭＶフィルタ処理は、この特性を利用して大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する。

大域的周波数分布を算出する方法として、上述した方法のほか、例えば、上述のＳＭＶフィルタ処理とＬ２ノルム正則化項を組み合わせた最小二乗処理により大域的周波数分布を算出する方法、三次元画像を低次多項式でフィッティングすることにより周波数分布を算出する方法、などがある。本実施形態では、こうした他の方法を用いてもよい。

［オフセット位相算出：Ｓ１２０６］
次に、オフセット位相算出部９０５において、低解像マルチエコー複素画像ｉ（ｎ）の各エコー画像から静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}による位相周りを除去することにより、各エコーのオフセット位相分布ｐ_{ｏｆｆｓｅｔ}（ｎ）を算出する。

ここで、低解像マルチエコー複素画像ｉ（ｎ）から静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}による位相周りを除去した複素画像をｉ’（ｎ）とすると、ｉ’（ｎ）は式（１５）を用いて式（２１）で表される。

したがって、各エコーのオフセット位相分布ｐ_{ｏｆｆｓｅｔ}（ｎ）は式（２２）より算出することができる。

式（２１）、（２２）に示すように、複素画像ｉ’（ｎ）の位相成分には静磁場不均一周波数ｆ_{ｉｎｈｏｍｏ}による位相周りが除去されていることが分かる。

［加重平均：Ｓ１２０７］
次に、加重平均部３３５において、式（２３）に示すように、各エコーのオフセット位相分布ｐ_{ｏｆｆｓｅｔ}（ｎ）を加重平均することで一つのオフセット位相分布ｐ_{ｏｆｆｓｅｔ}を算出する。

式中、ｗ（ｎ）は重みであり、前掲の式（９）や式（１０）に示す重みなど、任意の重みを用いてもよい。なお、この加重平均処理は式（２４）に示すように複素画像上で加重平均したのちに、位相成分を算出する方法でもよい。

この加重平均処理によりオフセット位相分布ｐ_{ｏｆｆｓｅｔ}のＳＮＲが改善するため、算出精度が向上する。

なお、図８に示す処理フローにおいて、Ｓ１２０５と、Ｓ１２０６およびＳ１２０７の処理の順番は問わない。すなわち、大域的周波数分布を算出した後（Ｓ１２０５）、加重平均し（Ｓ１２０７）、オフセット位相分布を算出してもよい（Ｓ１２０６）。

以上の位相分布分離処理において算出した低解像大域的周波数分布と低解像オフセット位相分布とをもとの解像度に変換した後（図６：Ｓ１１０３）、これらを計測した複素画像Ｉ（ｎ）から除去し（図６：Ｓ１１０４）、局所的周波数分布ｆ_{ｌｏｃａｌ}を算出する。さらに局所的周波数分布を用いて、磁場分布や、定量的磁化率分布、磁化率強調画像などの分布画像を取得することができる（図４：Ｓ１００４）。

本実施形態によれば、第一の実施形態と同様に、低解像度化した少ない画素数の分布を用いてオフセット位相分布と大域的周波数分布の分離を行うので演算時間が短縮され、さらに、複素画像から大域的周波数分布を除去する処理は、低解像度の段階で分離し高解像度に戻したオフセット位相分布と大域的周波数分布をもとの複素画像から除去するだけでよいので、さらに、演算時間が短縮される。また大域的周波数分布だけでなくオフセット位相分布も別に分離することで、背景除去処理（大域的周波数分布除去）の精度を高めることができる。さらに、エコー毎に得られる低解像オフセット位相分布や局所的周波数分布を加重加算することで、最終的に得られる局所的周波数分布のＳＮＲを向上することができる。

＜第三の実施形態＞
上述した第一の実施形態では、図８及び図９に示したように、位相分布分離部３３２は、まず低解像複素画像ｉ（ｎ）の位相差を求め、オフセット位相分布を除去してから位相アンラップ及び周波数変換を行って低解像大域的周波数分布を算出し、さらに低解像オフセット位相分布を算出した。

本実施形態では、位相分布分離部３３２は、信号フィッティング処理を用いることで、低解像複素画像ｉ（ｎ）から、一度に、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を分離する。

本実施形態においても、図６のフローチャートに示すように、低解像度化処理（Ｓ１１０１）が行われ、また位相分布分離処理（Ｓ１１０２）に続いて、高解像度化処理（Ｓ１１０３）、位相分布除去処理（Ｓ１１０４）等が行われることは第一及び第二の実施形態と同様である。

以下、本実施形態の、信号フィッティング処理を用いた位相分布分離処理（Ｓ１１０２）について図１０及び図１１を参照して説明する。図１０は位相分布分離部３３２の処理の概要を示す図、図１１は、位相分布分離部３３２の機能を示す機能ブロック図である。

位相分布分離部３３２は、図１０に示すように、低解像変換部３３１にて複素画像Ｉ（ｎ）を低解像度化した低解像複素画像ｉ（ｎ）において、ＧｒＥ法の信号モデルを用いたフィッティングにより、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}とオフセット位相分布ｐ_{ｏｆｆｓｅｔ}を算出する（Ｓ１３０１）。信号フィッティングで算出したオフセット位相成分をオフセット位相分布ｐ_{ｏｆｆｓｅｔ}として設定する（Ｓ１３０２）。そして、信号フィッティングで算出した静磁場不均一周波数を周波数アンラップして静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}を算出する（Ｓ１３０３）。大域的周波数分布算出処理（周波数分布分離処理）により、周波数アンラップした静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}から大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する（Ｓ１３０４）。以上の手順により、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を個別に処理することなくそれぞれ算出することができる。

これらを実現するため、本実施形態の位相分布分離部３３２は、図１１に示すように、低解像複素画像ｉ（ｎ）をＧｒＥ法の信号モデルに信号フィッティングさせる信号フィッティング部９１１と、信号フィッティング処理により得られたオフセット位相成分をオフセット位相分布ｐ_{ｏｆｆｓｅｔ}として設定するオフセット位相分布設定部９１２と、信号フィッティングで算出した静磁場不均一周波数分布を周波数アンラップして静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}を算出する周波数アンラップ部９１３と、周波数アンラップした静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}から大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する大域的周波数分布算出部９１４と、を備える。

以下、本変形例の位相分布分離部３３２の各部が行う処理の具体的な内容を、図１０を参照して詳述する。

［フィッティング：Ｓ１３０１］
信号フィッティング部９１１において、低解像度マルチエコー複素画像ｉ（ｎ）をＧｒＥ法の信号モデルに信号フィッティングさせる。ＧｒＥ法で計測した画素内の信号モデルＳ（ｎ）は、画素内のプロトン密度分布をＭ_０、見かけの横緩和速度分布をＲ_２ ^＊とすると、式（２５）で表される。

低解像マルチエコー複素画像ｉ（ｎ）を式（２５）の信号モデルＳ（ｎ）に最小二乗フィッティングさせることで、プロトン密度分布Ｍ_０、見かけの横緩和速度分布Ｒ_２ ^＊、一時的な静磁場不均一周波数分布ｆ’_{ｉｎｈｏｍｏ}、オフセット位相分布ｐ_{ｏｆｆｓｅｔ}をそれぞれ算出することができる。信号フィッティングには最小二乗法の他、最小二乗法に正則化項と呼ばれる制約項を付加した正則化項付き最小二乗法など、公知の方法を用いて信号フィッティングしてもよい。

［オフセット位相分布：Ｓ１３０２］
次に、Ｓ１３０２では、オフセット位相分布設定部９１２において、信号フィッティング処理により得られたオフセット位相成分をオフセット位相分布ｐ_{ｏｆｆｓｅｔ}として設定する。

［周波数アンラップ：Ｓ１３０３］
周波数アンラップ部９１３において、信号フィッティングで算出した一時的な静磁場不均一周波数分布ｆ’_{ｉｎｈｏｍｏ}に対し、周波数アンラップして静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}を算出する。

ここで、エコー間隔（エコー時間差分）ΔＴＥの逆数をｆ_ｎｙｑとすると、一時的な静磁場不均一周波数分布ｆ’_{ｉｎｈｏｍｏ}の一部の領域では、−ｆ_ｎｙｑ／２〜＋ｆ_ｎｙｑ／２の範囲を超えた周波数が−ｆ_ｎｙｑ／２〜＋ｆ_ｎｙｑ／２の範囲内に折り返される。これらの折り返しを補正する処理が周波数アンラップであり、本実施形態では、例えば領域拡大法などを用い、−ｆ_ｎｙｑ／２〜＋ｆ_ｎｙｑ／２の範囲に折り返されている周波数を補正する。このような周波数アンラップを行うことで撮像部位（例えば、頭部など）の全領域において正確な位相を求めることができる。

［大域的周波数分布算出：Ｓ１３０４］
次に、大域的周波数分布算出部９１４において、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}から、生体形状などに起因して生じる局所的周波数変化と大域的周波数変化とを分離し、大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する。大域的周波数分布算出部９１４が行う処理は、第二の実施形態の大域的周波数分布算出部９０４の処理（図７：Ｓ１２０５）と同様であり、例えば、ＳＭＶフィルタ処理や、その他の手法が採用できる。

以上の手順により、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を一度に算出することができる。その後、算出した低解像大域的周波数分布と低解像オフセット位相分布とを計測した複素画像と同じ解像度に変換した後（図６：Ｓ１１０３）、これらを計測した複素画像から除去し（図６：Ｓ１１０４）、局所的周波数分布局所的周波数分布ｆ_{ｌｏｃａｌ}を算出する。さらに局所的周波数分布を用いて、磁場分布や、定量的磁化率分布、磁化率強調画像などの分布画像を取得することができる（図４：Ｓ１００４）。

本実施形態によれば、低解像度化したデータでフィッティングするため、従来のフィッティング法に比べて演算時間を短縮できる。また、低解像画像であるため、ノイズが低減しており、信号フィッティングの精度が高まる。また、本実施形態によれば、演算時間のかかる大域的周波数分布算出処理が一回で済むため、従来の個別処理法に比べて演算時間を短縮できる。また、オフセット位相分布ｐ_{ｏｆｆｓｅｔ}分離した上で大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出しているため、大域的周波数分布ｆ_{ｇｌｏｂａｌ}の算出精度が向上する。

＜第四の実施形態＞
次に、本発明の第四の実施形態について説明する。第三の実施形態では、対象とする物質が一つの場合（例えば水信号のみ）に関しての実施形態である。一方、本実施形態は、対象とする部位に、第一の物質（例えば水信号）と第一の物質と共鳴周波数が異なる第二の物質（例えば脂肪信号）が混在する場合について、少なくとも２つ以上の異なるＴＥのマルチエコー複素画像から局所的周波数分布を算出する。

本実施形態で用いるＭＲＩ装置１００は、基本的に第一の実施形態と同様の構成を有する。また図４及び図６に示す処理の概要も共通であるが、位相分布分離処理（図６：Ｓ１１０２）の内容が異なる。以下、本実施形態について、第一の実施形態と異なる位相分離部３３２の構成に主眼をおいて説明する。

本実施形態の信号フィッティング処理を用いた位相分布分離処理（Ｓ１１０２）について図１２を用いて説明する。なお本実施形態では第一の物質を水、第二の物質を脂肪として説明する。

位相分布分離部３３２は、低解像変換部３３１にて複素画像Ｉ（ｎ）を低解像度化した低解像複素画像ｉ（ｎ）において、ＧｒＥ法の信号モデルを用いたフィッティングにより、水の信号強度Ｗと脂肪の信号強度Ｆ、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}とゼロ次位相分布を算出する（Ｓ１４０１）。そして、信号フィッティングで算出したゼロ次位相成分をオフセット位相分布ｐ_{ｏｆｆｓｅｔ}として設定する（Ｓ１４０２）。オフセット位相分布ｐ_{ｏｆｆｓｅｔ}と、水と脂肪の周波数差Ｆ_ｆａｔに起因してＴＥ_ｎによって変化する位相成分p_ｆａｔ（ｎ）を合成して合成オフセット位相分布ｐ_ａｄｄ（ｎ）を算出する（Ｓ１４０３）。そして、信号フィッティングで算出した静磁場不均一周波数を周波数アンラップして静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}を算出する（Ｓ１４０４）。大域的周波数算出処理により、周波数アンラップした静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}から大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する（Ｓ１４０５）。以上の手順により、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と、低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}（水と脂肪の周波数差に起因する位相差p_ｆａｔを含む）をそれぞれ算出することができる。

これらを実現するため、本実施形態の位相分布分離部３３２は、図１３に示すように、低解像複素画像ｉ（ｎ）を、水と脂肪を含んだＧｒＥ法の信号モデルに信号フィッティングさせる信号フィッティング部９２１と、信号フィッティング処理により得られたオフセット位相成分をオフセット位相分布ｐ_{ｏｆｆｓｅｔ}として設定するオフセット位相分布設定部９２２と、オフセット位相分布ｐ_{ｏｆｆｓｅｔ}と、水と脂肪の周波数差f_ｆａｔに起因する位相成分p_ｆａｔ（ｎ）を合成して、合成オフセット位相分布ｐ_ａｄｄ（ｎ）を算出する合成オフセット位相分布算出部９２３と、信号フィッティングで算出した静磁場不均一周波数分布を周波数アンラップして静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}を算出する周波数アンラップ部９２４と、大域的周波数分布算出処理により、周波数アンラップした静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}から大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する大域的周波数分布算出部９２５と、を備える。

以下、図１２の各処理について説明する。
［フィッティング：Ｓ１４０１］
まずＳ１４０１では、信号フィッティング部９２１において、低解像複素画像ｉ（ｎ）をＧｒＥ法の信号モデルを用いたフィッティングにさせて、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}とオフセット位相分布ｐ_{ｏｆｆｓｅｔ}を算出する。ここで、画素内の水の信号強度をＷ、脂肪の信号強度をＦ、水と脂肪の共鳴周波数差をｆ_ｆａｔとすると、ＧｒＥ法で計測した画素内の信号モデルＳ（ｎ）は式（２６）で表される。

計測した低解像マルチエコー複素画像ｉ（ｎ）を式（２６）の信号モデルＳ（ｎ）に最小二乗フィッティングさせることで、画素内の水の信号強度Ｗ、脂肪の信号強度Ｆ、一時的な静磁場不均一周波数分布ｆ’_{ｉｎｈｏｍｏ}、オフセット位相分布ｐ_{ｏｆｆｓｅｔ}をそれぞれ算出することができる。信号フィッティングには最小二乗法の他、最小二乗法に正則化項と呼ばれる制約項を付加した正則化項付き最小二乗法など、公知の方法を用いて信号フィッティングしてもよい。また式（２６）の｛Ｗ＋Ｆ（ｊ２πｆ_ｆａｔＴＥ_ｎ）｝は水と脂肪の共鳴周波数差ｆ_ｆａｔに起因する位相成分p_ｆａｔ（ｎ）であり、ＷとＦを用いて算出することができる。

［オフセット位相分布設定：Ｓ１４０２］
オフセット位相分布設定部９２２において、信号フィッティング処理により得られた位相成分をオフセット位相分布ｐ_{ｏｆｆｓｅｔ}として設定する。

［合成オフセット位相分布算出：Ｓ１４０３］
また、水と脂肪の共鳴周波数差に起因する位相成分p_ｆａｔ（ｎ）は、オフセット位相分布と同様に局所的周波数分布ｆ_{ｌｏｃａｌ}を算出するために除去すべき成分である。そこでＳ１４０３では、合成オフセット位相分布算出部９２３において、まず式（２７）に示すように、オフセット位相分布ｐ_{ｏｆｆｓｅｔ}と、水と脂肪の周波数差f_ｆａｔに起因する位相成分p_ｆａｔ（ｎ）を合成した位相分布（合成オフセット位相分布という）ｐ_ａｄｄ（ｎ）を算出する。

［周波数アンラップ：Ｓ１４０４］
周波数アンラップ部９２４において、信号フィッティングで算出した一時的な静磁場不均一周波数分布ｆ’_{ｉｎｈｏｍｏ}に対し、周波数アンラップして静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}を算出する。周波数アンラップの手法は、第三の実施形態と同様であり、−ｆ_ｎｙｑ／２〜＋ｆ_ｎｙｑ／２の範囲を超えて、−ｆ_ｎｙｑ／２〜＋ｆ_ｎｙｑ／２の範囲内に折り返された周波数成分を補正する。なお「ｆ_ｎｙｑ」は、エコー間隔（エコー時間差分）ΔＴＥの逆数である。

［周波数分布分離：Ｓ１４０５］
次に、大域的周波数分布算出部（周波数分布分離部）９２５において、静磁場不均一周波数分布ｆ_{ｉｎｈｏｍｏ}から、生体形状などに起因して生じる局所的周波数変化と大域的周波数変化とを分離し、大域的周波数分布ｆ_{ｇｌｏｂａｌ}を算出する。この算出手法は、第二の実施形態の周波数分布分離部９０４の処理（図８：Ｓ１２０５）と同様であり、例えば、ＳＭＶフィルタ処理や、その他の手法が採用できる。

以上の手順により、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と、オフセット位相分布ｐ_{ｏｆｆｓｅｔ}と脂肪による位相変化ｐ_ｆａｔを含む低解像合成オフセット位相分布ｐ_ａｄｄを算出することができる。その後、高解像度変換部３３３にて、低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像合成オフセット位相分布ｐ_ａｄｄを高解像度化したのち、位相除去部３３４にて、各エコーの局所的周波数分布を算出し、加重平均部３３５にて、最終的な局所周波数分布を算出する。

その後、図４に示したように、図１３に示す位相分布分離処理において算出した低解像大域的周波数分布と低解像合成オフセット位相分布とをもとの解像度に変換した後（図６：Ｓ１１０３）、これらを計測した複素画像Ｉ（ｎ）から除去し（図６：Ｓ１１０４）、局所的周波数分布ｆ_{ｌｏｃａｌ}を算出する。さらに局所的周波数分布を用いて、磁場分布や、定量的磁化率分布、磁化率強調画像などの分布画像を取得することができる（図４：Ｓ１００４）。

本実施形態によれば、水と脂肪が混在するような部位においても短時間かつ高精度に低解像大域的周波数分布ｆ_{ｇｌｏｂａｌ}と低解像オフセット位相分布ｐ_{ｏｆｆｓｅｔ}を算出することができる。また最終的に局所周波数分布を求める際に不要となる水と脂肪との共鳴周波数差による位相成分ｐ_ｆａｔも算出することができ、高精度で局所周波数分布を算出することができる。その結果、水と脂肪が混在する部位においても，脂肪による位相周りの影響を除去したＱＳＭやＳＷＩが算出できる。

１００：ＭＲＩ装置、１０１：被検体、１０２：静磁場コイル、１０３：傾斜磁場コイル、１０４：シムコイル、１０５：送信コイル、１０６：受信コイル、１０７：送信機、１０８：受信機、１０９：計算機、１１０：ディスプレイ、１１１：外部記憶装置、１１２：傾斜磁場用電源部、１１３：シム用電源部、１１４：シーケンス制御装置、１１５：入力装置、１２０：ＭＲＩ装置、１３０：ＭＲＩ装置、３１０：計測制御部、３２０：画像再構成部、３３０：局所的周波数分布算出部、３３１：低解像変換部（第一の解像度変換部）、３３２：位相分布分離部、３３３：高解像変換部（第二の解像度変換部）、３３４：位相除去部、３３５：加重平均部、３４０：分布画像算出部、９０１：位相差算出部、９０２：位相アンラップ部、９０３：周波数変換部、９０４：大域的周波数分布算出部（周波数分布分離部）、９０５：オフセット位相算出部、９１１：信号フィッティング部、９１２：オフセット位相分布設定部、９１３：周波数アンラップ部、９１４：大域的周波数分布算出部（周波数分布分離部）、９２１：信号フィッティング部、９２２：オフセット位相分布設定部、９２３：合成オフセット位相分布算出部、９２４：周波数アンラップ部、９２５：大域的周波数分布算出部（周波数分布分離部）

Claims (12)

1. 静磁場内に配置された被検体に高周波磁場パルスを送信する送信部と、
   前記被検体が発生する核磁気共鳴信号を受信する受信部と、
   静磁場に傾斜磁場を与える傾斜磁場発生部と、
   前記受信した核磁気共鳴信号に演算を施す計算機と、を備えた磁気共鳴イメージング装置であって、
   前記計算機は、
   複数の異なるエコー時間で取得した核磁気共鳴信号から複数の複素画像を生成する画像再構成部と、
   前記複数の複素画像を、それぞれ、当該複素画像より低解像度に変換する第一の解像度変換部と、
   前記第一の解像度変換部で処理された低解像度画像から、大域的周波数分布とオフセット位相分布とを分離する位相分布分離部と、
   前記位相分布分離部によって分離された大域的周波数分布及びオフセット位相分布を、前記複数の複素画像と同じ解像度に変換する第二の解像度変換部と、
   前記第二の解像度変換部で処理された大域的周波数分布及びオフセット位相分布と、前記複数の複素画像とを用いて、局所的周波数分布を算出する局所的周波数分布算出部とを備えたことを特徴とする磁気共鳴イメージング装置。
2. 請求項１に記載の磁気共鳴イメージング装置であって、
   前記位相分布分離部は、
   複数の前記低解像度画像について、画像間の位相差を算出する位相差算出部と、
   前記位相差算出部が算出した位相差の位相折り返しを除去する折り返し除去部と、
   前記折り返し除去部により折り返しが除去された位相差を、周波数に換算して周波数分布を算出する周波数変換部と、
   前記周波数分布算出部が算出した周波数分布を、大域的周波数分布と局所的周波数分布とに分離する周波数分布分離部と、
   前記低解像度画像と前記周波数変換部が算出した周波数分布とを用いて、オフセット位相分布を算出するオフセット位相分布算出部と、
   を備えることを特徴とする磁気共鳴イメージング装置。
3. 請求項１に記載の磁気共鳴イメージング装置であって、
   前記位相分布分離部は、
   複数の前記低解像度画像を、それぞれ、計測時の信号論理式にフィッティングさせて周波数分布とオフセット位相分布とを算出するフィッティング部と、
   前記周波数分布の折り返しを除去する折り返し除去部と、
   前記折り返し除去部により折り返しが除去された周波数分布を、大域的周波数分布と局所的周波数分布とに分離する周波数分布分離部と、
   を備えることを特徴とする磁気共鳴イメージング装置。
4. 請求項１に記載の磁気共鳴イメージング装置であって、
   前記位相分布分離部は、
   第一の物質の信号と、当該第一の物質とは共鳴周波数が異なる第二の物質の信号とを含む信号論理式を設定し、複数の前記低解像度画像を、それぞれ、前記信号論理式にフィッティングさせて、前記第一の物質と前記第二の物質との共鳴周波数差に起因する位相分布と、周波数分布と、オフセット位相分布とを算出するフィッティング部と、
   前記フィッティング部により算出された周波数分布を、大域的周波数分布と局所的周波数分布とに分離する周波数分布分離部と、
   前記共鳴周波数差に起因する位相分布と前記オフセット位相分布とを合成し、合成オフセット位相分布を算出する合成オフセット位相分布算出部とを備え、
   前記第二の解像度変換部は、前記大域的周波数分布と前記合成オフセット位相分布を、前記複数の複素画像と同じ解像度に変換し、前記局所的周波数分布算出部は、前記第二の解像度変換部で処理された大域的周波数分布及び合成オフセット位相分布と複数の複素画像とを用いて、局所的周波数分布を算出することを特徴とする磁気共鳴イメージング装置。
5. 請求項４に記載の磁気共鳴イメージング装置であって、
   前記第一の物質及び前記第二の物質が、水及び脂肪であることを特徴とする磁気共鳴イメージング装置。
6. 請求項１に記載の磁気共鳴イメージング装置であって、
   前記計算機は、複数の前記複素画像毎に算出した前記大域的周波数分布及び／又は局所的周波数分布を加重平均する加重平均部をさらに備えることを特徴とする磁気共鳴イメージング装置。
7. 請求項１に記載の磁気共鳴イメージング装置であって、
   前記第二の解像度変換部は、大域的周波数分布及び／又はオフセット位相分布を平滑化する平滑部を備えることを特徴とする磁気共鳴イメージング装置。
8. 請求項１に記載の磁気共鳴イメージング装置であって、
   前記局所的周波数分布算出部が算出した局所的周波数分布を用いて、局所磁場分布を算出する局所磁場算出部と、
   前記局所磁場算出部が算出した局所磁場と、磁場と磁化率との関係式とを用いて磁化率分布を算出する磁化率分布算出部と、を
   さらに備えることを特徴とする磁気共鳴イメージング装置。
9. 請求項１に記載の磁気共鳴イメージング装置であって、
   前記局所的周波数分布算出部が算出した局所的周波数分布を用いて、磁化率強調マスクを作成するマスク作成部と、
   前記複数の複素画像の絶対値成分に、前記マスク生成部が作成した磁化率強調マスクを乗じて磁化率強調画像を作成する磁化率強調画像作成部と、を
   さらに備えることを特徴とする磁気共鳴イメージング装置。
10. 請求項１に記載の磁気共鳴イメージング装置であって、
    前記計算機を、独立した画像処理装置として有することを特徴とする磁気共鳴イメージング装置。
11. 磁気共鳴イメージング装置により異なるエコー時間で取得した複数の複素画像を用いて、局所周波数分布を算出する画像処理方法であって、
    前記複数の複素画像を、それぞれ、当該複素画像より低解像度の画像データに変換し、
    低解像度の複素画像から、大域的周波数分布とオフセット位相分布とを分離し、
    前記大域的周波数分布及びオフセット位相分布を前記複数の画像データと同じ高解像度のデータに変換し、
    高解像度に変換された大域的周波数分布及びオフセット位相分布と前記複数の複素画像とを用いて、局所的周波数分布を算出することを特徴とする画像処理方法。
12. 磁気共鳴イメージング装置により異なるエコー時間で取得した複数の複素画像を用いて、局所周波数分布を算出する画像処理方法であって、
    前記複数の複素画像を、それぞれ、当該複素画像より低解像度の画像データに変換し、
    第一の物質の信号と、当該第一の物質とは共鳴周波数が異なる第二の物質の信号とを含む信号論理式を設定し、複数の前記低解像度画像を、それぞれ、前記信号論理式にフィッティングさせて、前記第一の物質と前記第二の物質との共鳴周波数差に起因する位相分布と、周波数分布と、オフセット位相分布とを算出し、
    前記周波数分布を、大域的周波数分布と局所的周波数分布とに分離し、
    前記共鳴周波数差に起因する位相分布と前記オフセット位相分布とを合成し、合成位相分布を算出し、
    前記周波数分布から分離された大域的周波数分布と前記合成位相分布とを、前記複数の複素画像と同じ解像度に変換し、
    前記複数の複素画像と同じ解像度に変換された大域的周波数分布及び合成位相分布と前記複数の複素画像とを用いて、局所的周波数分布を算出することを特徴とする画像処理方法。
    
    

Images