(57)【要約】  ＭＲＩを用いて磁化率分布を算出するにあたり、簡易な処理で、アーチファクトおよびノイズを低減し、算出する磁化率分布の精度を向上させる技術を提供する。この技術では、まず、ＭＲＩで取得した複素画像から位相画像を介して磁場分布を得る。次いで磁場と磁化率との関係式に基づく制約条件下で、磁場分布から算出した磁化率分布に対して平滑化処理を行うことを繰り返す。具体的には、予め定めた初期磁化率分布と磁場分布とにより定義した誤差関数を繰り返し演算により最小化する最小化処理により前記磁化率分布を算出する。このとき、最小化処理中に前記誤差関数に基づいて算出される更新磁化率分布を、当該更新磁化率分布が算出される毎に平滑化する。
2016076076.tif 000002
【特許請求の範囲】
【請求項１】
高周波磁場パルスの照射に応じて被検体から発生する核磁気共鳴信号を複素信号として計測する計測制御部と、
前記複素信号から画素値が複素数である複素画像を再構成する画像再構成部と、
前記複素画像から位相画像を生成し、当該位相画像から磁場分布を算出する磁場分布算出部と、
前記磁場分布から磁化率分布を算出する磁化率分布算出部と、を備え、
前記磁化率分布算出部は、予め定めた初期磁化率分布と前記磁場分布とを用いて定義した誤差関数を繰り返し演算により最小化する最小化処理を行って前記磁化率分布を算出するとともに、当該最小化処理中に前記誤差関数に基づいて算出される更新磁化率分布を、当該更新磁化率分布が算出される毎に平滑化すること
を特徴とする磁気共鳴イメージング装置。
【請求項２】
請求項１記載の磁気共鳴イメージング装置であって、
前記磁化率分布算出部は、
前記誤差関数を用いて前記初期磁化率分布から前記更新磁化率分布を算出する更新磁化率分布算出部と、
前記更新磁化率分布を平滑化して合成磁化率分布を生成する平滑化処理部と、
前記繰返し演算の終了判定を行うとともに、前記初期磁化率分布を前記合成磁化率分布に更新する判定部と、を備え、
前記平滑化処理部は、
前記更新磁化率分布に予め定めた平滑化フィルタを適用し、処理後磁化率分布を生成するフィルタ処理部と、
前記更新磁化率分布と前記処理後磁化率分布とを合成し、合成磁化率分布を生成する合成部と、を備え、
前記判定部は、終了と判定した際の前記合成磁化率分布を、前記磁化率分布とすること
を特徴とする磁気共鳴イメージング装置。
【請求項３】
請求項１または２記載の磁気共鳴イメージング装置であって、
前記誤差関数は、前記磁場分布と前記初期磁化率分布に畳み込み演算を行った結果とにより定まる関数に重みを付与して設定され、
前記重みは、前記複素画像から得られる位相画像に基づいて算出されること
を特徴とする磁気共鳴イメージング装置。
【請求項４】
請求項２記載の磁気共鳴イメージング装置であって、
前記更新磁化率分布算出部は、前記誤差関数の勾配を算出し、算出した勾配を用いて、前記更新磁化率分布を算出すること
を特徴とする磁気共鳴イメージング装置。
【請求項５】
請求項２記載の磁気共鳴イメージング装置であって、
前記合成部は、前記更新磁化率分布および前記処理後磁化率分布をそれぞれｋ空間データに変換し、前記変換後のｋ空間データそれぞれに重みを付けて加算することにより、前記合成磁化率分布を生成すること
を特徴とする磁気共鳴イメージング装置。
【請求項６】
請求項５記載の磁気共鳴イメージング装置であって、
前記更新磁化率分布のｋ空間データにおいて、マジックアングル領域近傍のデータに付与する重みは、他の領域のデータに付与する重みより小さく設定され、
前記処理後磁化率分布のｋ空間データにおいて、マジックアングル領域近傍のデータに付与する重みは、他の領域のデータに付与する重みより大きく設定されること
を特徴とする磁気共鳴イメージング装置。
【請求項７】
請求項５記載の磁気共鳴イメージング装置であって、
前記重みは、１つの前記ｋ空間データ内の全データに一様に設定されること
を特徴とする磁気共鳴イメージング装置。
【請求項８】
請求項２記載の磁気共鳴イメージング装置であって、
前記フィルタは、エッジ保存平滑化フィルタであること
を特徴とする磁気共鳴イメージング装置。
【請求項９】
請求項２記載の磁気共鳴イメージング装置であって、
前記初期磁化率分布は、全ての値が０であること
を特徴とする磁気共鳴イメージング装置。
【請求項１０】
請求項２記載の磁気共鳴イメージング装置であって、
前記判定部は、前記更新磁化率分布が算出されると、当該更新磁化率分布と前記初期磁化率分布とを用いて前記終了判定を行うこと
を特徴とする磁気共鳴イメージング装置。
【請求項１１】
請求項２記載の磁気共鳴イメージング装置であって、
前記判定部は、前記合成磁化率分布が算出されると、当該合成磁化率分布と前記初期磁化率分布とを用いて前記終了判定を行うこと
を特徴とする磁気共鳴イメージング装置。
【請求項１２】
高周波磁場パルスの照射に応じて被検体から発生する核磁気共鳴信号を収集し、再構成することにより得た複素画像から位相画像を得、
前記位相画像から磁場分布を算出し、
予め定めた初期磁化率分布と前記磁場分布とにより定義した誤差関数を、繰り返し演算により最小化する最小化処理を行い、磁化率分布を算出し、
前記最小化処理中に、
前記誤差関数から算出される、前記初期磁化率分布を更新した更新磁化率分布に対して平滑化を行い、合成磁化率分布を得、
前記初期磁化率分布を当該合成磁場分布に更新すること
を特徴とする磁気共鳴イメージング装置における定量的磁化率マッピング方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、磁気共鳴イメージング（ＭＲＩ：Ｍａｇｎｅｔｉｃ  Ｒｅｓｏｎａｎｃｅ  Ｉｍａｇｉｎｇ）技術に関する。特に、取得した画像を用いて生体内の磁化率分布を算出する画像処理技術に関する。
【背景技術】
【０００２】
磁気共鳴イメージング装置は、静磁場内に置かれた水素原子核（プロトン）が特定の周波数の高周波磁場に共鳴する核磁気共鳴現象を利用した非侵襲な医用画像診断装置である。核磁気共鳴信号はプロトン密度や緩和時間など様々な物性値によって変化するため、ＭＲＩで得られる画像は生体組織の構造や組成、細胞性状などの様々な生体情報を描出できる。【０００３】
近年、ＭＲＩで測定できる物性値として、生体組織間の磁化率差が注目されている。磁化率とは、静磁場中の物質が磁気分極（磁化）する度合いを表した物性値である。生体内には、静脈血中のデオキシヘモグロビンや鉄たんぱく質などの常磁性体と、生体組織の大部分を占める水や石灰化の基となるカルシウムなどの反磁性体とがある。生体組織間の磁化率の差である磁化率差を定量的に画像化することで、脳虚血疾患の診断やがんの放射線治療効果の予測、神経変性疾患の鑑別に適用できる可能性がある。【０００４】
ＭＲＩを用いて生体組織間の磁化率差を画像化する方法は、定量的磁化率マッピング（ＱＳＭ：Ｑｕａｎｔｉｔａｔｉｖｅ  Ｓｕｓｃｅｐｔｉｂｉｌｉｔｙ  Ｍａｐｐｉｎｇ）法と呼ばれる。ＱＳＭ法は、計測したＭＲ画像の位相情報から生体組織間の磁化率差によって生じる空間的な磁場変化を算出し、磁場と磁化率の関係式を用いて磁化率分布を算出する手法である。【０００５】
ところが、磁場分布は、磁化率分布を空間的に畳み込み積分したものであるため、磁場分布から磁化率分布を求める場合、畳み込み積分の逆問題となり、一意に解が得られない。【０００６】
一般に、最小二乗法によって、磁場分布から磁化率分布を求める。このとき、誤差関数を導入し、これを最小化する値を解とする。代表的な方法には、例えば、磁場分布および点ダイポール磁場のｋ空間上の演算により磁化率分布を算出するＴＫＤ（Ｔｒｕｎｃａｔｅｄ−ｂａｓｅｄ  Ｋ−ｓｐａｃｅ  Ｄｉｖｉｓｉｏｎ）法（例えば、非特許文献１参照）、ＴＫＤ法で算出した磁化率分布と、閾値処理により微細構造を抽出した磁化率分布とを繰り返し演算により合成するＩｔｅｒａｔｉｖｅ  ＳＷＩＭ（Ｓｕｓｃｅｐｔｉｂｉｌｉｔｙ  Ｗｅｉｇｈｔｅｄ  Ｉｍａｇｉｎｇ  ａｎｄ  Ｍａｐｐｉｎｇ）法（例えば、非特許文献２、特許文献１参照）、正則化付き最小二乗法を用いたＭＥＤＩ法（Ｍｏｒｐｈｏｌｏｇｙ  ｅｎａｂｌｅｄ  ｄｉｐｏｌｅ  ｉｎｖｅｒｓｉｏｎ）法（例えば、非特許文献３、特許文献２参照）がある。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】米国特許出願公開２０１１／０２６２０１７号明細書
【特許文献２】米国特許出願公開２０１１／００４４５２４号明細書
【非特許文献】
【０００８】
【非特許文献１】Ｓｈｍｕｅｌｉ  Ｋ他、  ”Ｍａｇｎｅｔｉｃ  Ｓｕｓｃｅｐｔｉｂｉｌｉｔｙ  Ｍａｐｐｉｎｇ  ｏｆ  Ｂｒａｉｎ  Ｔｉｓｓｕｅ  Ｉｎ  Ｖｉｖｏ  Ｕｓｉｎｇ  ＭＲＩ  Ｐｈａｓｅ  Ｄａｔａ”、  Ｍａｇｎｅｔｉｃ  Ｒｅｓｏｎａｎｃｅ  ｉｎ  Ｍｅｄｉｃｉｎｅ、２００９年、６２巻、１５１０−１５２２頁
【非特許文献２】Ｔａｎｇ  Ｊ他、  ”Ｉｍｐｒｏｖｉｎｇ  Ｓｕｓｃｅｐｔｉｂｉｌｉｔｙ  Ｍａｐｐｉｎｇ  Ｕｓｉｎｇ  ａ  Ｔｈｒｅｓｈｏｌｄ−Ｂａｓｅｄ  Ｋ−Ｓｐａｃｅ／Ｉｍａｇｅ  Ｄｏｍａｉｎ  Ｉｔｅｒａｔｉｖｅ  Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ  Ａｐｐｒｏａｃｈ”、  Ｍａｇｎｅｔｉｃ  Ｒｅｓｏｎａｎｃｅ  ｉｎ  Ｍｅｄｉｃｉｎｅ、２０１３年、６９巻、１３９６−１４０７頁
【非特許文献３】Ｌｉｕ  Ｔ他、”Ｍｏｒｐｈｏｌｏｇｙ  ｅｎａｂｌｅｄ  ｄｉｐｏｌｅ  ｉｎｖｅｒｓｉｏｎ  （ＭＥＤＩ）  ｆｒｏｍ  ａ  ｓｉｎｇｌｅ−ａｎｇｌｅ  ａｃｑｕｉｓｉｔｉｏｎ：  ｃｏｍｐａｒｉｓｏｎ  ｗｉｔｈ  ＣＯＳＭＯＳ  ｉｎ  ｈｕｍａｎ  ｂｒａｉｎ  ｉｍａｇｉｎｇ”、Ｍａｇｎｅｔｉｃ  Ｒｅｓｏｎａｎｃｅ  ｉｎ  Ｍｅｄｉｃｉｎｅ、２０１１年、６６巻、７７７−７８３頁
【発明の概要】
【発明が解決しようとする課題】
【０００９】
ＴＫＤ法は、逆計算が不能な領域を、０または任意の値に打ち切る。このため、ＴＫＤ法は、簡便かつ計算時間が短いという利点があるものの、マジックアングル領域近傍を打ち切ることによりアーチファクトや算出磁化率の確度の低下が生じるという課題がある。【００１０】
Ｉｔｅｒａｔｉｖｅ  ＳＷＩＭ法は、ＴＫＤ法で得た磁化率分布に対し、微細構造を維持可能なマスク処理し、磁化率分布を更新する処理を繰り返す。微細構造を維持した磁化率分布を算出できるものの、アーチファクトを完全には除去できないという課題がある。また、微細構造を抽出するためには、さまざまな閾値処理が必要なため、被検体ごとに適切な閾値を設定しないと微細構造を適切に抽出できないという課題がある。【００１１】
ＭＥＤＩ法は、磁場と磁化率との関係式を制約条件として、正則化項と呼ばれる平滑効果を与える項を繰り返し演算により最小化することで磁化率分布を算出する。このような正則化付き最小二乗法は、正則化パラメータの値に対応した大きさの信号値を打ち切る処理であるため、磁場分布に混入したノイズに応じて、適切に正則化項を設定する必要がある。しかしながら、通常、磁場分布に混入したノイズの情報を事前に得ることができないため、正則化パラメータを複数回変化させて算出した残差の変化曲線から適切な正則化パラメータを算出する。従って、ＭＥＤＩ法では、適切な正則化パラメータの設定のために膨大な計算時間が必要となる。【００１２】
本発明は、上記事情に鑑みて成されたもので、ＭＲＩを用いて磁化率分布を算出するにあたり、簡易な処理で、アーチファクトおよびノイズを低減し、算出する磁化率分布の確度および精度を向上させる技術を提供することを目的とする。
【課題を解決するための手段】
【００１３】
本発明は、ＭＲＩで取得した複素画像から位相画像を介して磁場分布を得る。磁場と磁化率との関係式に基づく制約条件下で、磁場分布から算出した磁化率分布に対して平滑化処理を行うことを繰り返す。具体的には、予め定めた初期磁化率分布と磁場分布とにより定義した誤差関数を繰り返し演算により最小化する最小化処理により前記磁化率分布を算出する。このとき、最小化処理中に前記誤差関数に基づいて算出される更新磁化率分布を、当該更新磁化率分布が算出される毎に平滑化する。
【発明の効果】
【００１４】
ＭＲＩを用いて磁化率分布を算出するにあたり、簡易な処理で、アーチファクトおよびノイズを低減し、算出する磁化率分布の確度および精度が向上する。
【図面の簡単な説明】
【００１５】
【図１】（ａ）は、本発明の実施形態の、垂直磁場方式の磁気共鳴イメージング装置の外観図、（ｂ）は、本発明の実施形態の、水平磁場方式の磁気共鳴イメージング装置の外観図、（ｃ）は、本発明の実施形態の、開放感を高めた磁気共鳴イメージング装置の外観図である。
【図２】本発明の実施形態におけるＭＲＩ装置の概略構成を示すブロック図である。
【図３】本発明の実施形態の計算機の機能ブロック図である。
【図４】本発明の実施形態のＱＳＭ処理のフローチャートである。
【図５】ＲＳＳＧ（ＲＦ−ｓｐｏｉｌｅｄ−Ｓｔｅａｄｙ−ｓｔａｔｅ  Ａｃｑｕｉｓｉｔｉｏｎ  ｗｉｔｈ  Ｒｅｗｏｕｎｄ  Ｇｒａｄｉｅｎｔ−Ｅｃｈｏ）シーケンスのパルスシーケンスを説明するための説明図である。
【図６】本発明の実施形態の磁場分布算出処理のフローチャートである。
【図７】本発明の実施形態の磁化率分布算出処理のフローチャートである。
【図８】本発明の実施形態の磁化率分布算出処理の流れを説明するための説明図である。
【図９】本発明の実施形態の変形例の磁化率分布算出処理のフローチャートである。
【図１０】本発明の実施形態の変形例の磁化率分布算出処理の流れを説明するための説明図である。
【発明を実施するための形態】
【００１６】
本発明を適用する実施形態の一例を説明する。以下、実施形態を説明するための全図において、特に断らない限り、同一機能を有するものは同一符号を付し、その繰り返しの説明は省略する。また、以下の記述により本発明が限定されるものではない。
【００１７】
［磁化率分布算出の原理］
本実施形態では、ＭＲＩ装置を用いて、予め定めたエコー時間の複素画像Ｉを撮影し、複素画像Ｉから生体組織間の磁化率差によって生じる磁場変化をとらえた磁場分布δを算出し、磁場分布δから、磁化率分布χを算出する。いわゆるＱＳＭ法により磁化率分布を算出する。これを実現するＭＲＩ装置の説明に先立ち、ＱＳＭ法において、複素画像Ｉから磁化率分布χを算出する手法について説明する。
【００１８】
ＱＳＭ法では、生体組織間の磁化率差によって生じる局所的な磁場変化を、Ｇｒａｄｉｅｎｔ  Ｅｃｈｏ（ＧｒＥ）法で撮像した位相画像から算出する。
【００１９】
位置ベクトルをｒとすると組織間の磁化率差によって生じる相対的な磁場変化（磁場分布）δ（ｒ）は、以下の式（１）で表される。
【数１】
2016076076.tif 000003
ここで、Ｐ（ｒ）は、位相画像、γは、プロトンの核磁気回転比、Ｂ0は、静磁場強度、ＴＥは、エコー時間をそれぞれ表す。
【００２０】
また、磁場分布δ（ｒ）は、静磁場に関するマクスウェル方程式より、生体内の磁化率分布χ（ｒ）を用いて以下の式（２）で表される。
【数２】
2016076076.tif 000004
ここで、αは、ベクトル（ｒ’−ｒ）と静磁場方向とのなす角度、ｄ（ｒ）は、点ダイポール磁場をそれぞれ表す。
【００２１】
式（２）に示すように、磁場分布δ（ｒ）は、磁化率分布χ（ｒ）と点ダイポール磁場ｄ（ｒ）との畳み込み積分で表される。したがって、式（２）の両辺をフーリエ変換することにより、式（２）は以下の式（３）に変換される。
【数３】
2016076076.tif 000005
ここで、ｋ＝（ｋx、ｋy、ｋz）は、ｋ空間上の位置ベクトル、Δ（ｋ）、Ｘ（ｋ）、Ｄ（ｋ）、は、磁場分布δ（ｒ）、磁化率分布χ（ｒ）、点ダイポール磁場ｄ（ｒ）のフーリエ成分をそれぞれ表す。
【００２２】
式（３）に示すように、磁化率分布のフーリエ成分Ｘ（ｋ）は、磁場分布のフーリエ成分Δ（ｋ）を点ダイポール磁場のフーリエ成分Ｄ（ｋ）で除算することによって算出できる。しかしながら、式（３）は、Ｄ（ｋ）＝０近傍の領域において、その逆数が発散してしまうため、直接的にＸ（ｋ）を算出することができない。
【００２３】
このＤ（ｋ）＝０となる領域はマジックアングルと呼ばれ、磁場方向に対しておよそ５４．７°の２倍の頂角を持つ逆双円錐領域となる。マジックアングルの存在により磁場分布から磁化率分布を算出するＱＳＭ法は、不良条件逆問題（ｉｌｌ−ｃｏｎｄｉｔｉｏｎｅｄ  ｉｎｖｅｒｓｅ  ｐｒｏｂｌｅｍ）に帰着される。
【００２４】
この問題を解く代表的な方法として、上述のように、ＴＫＤ法と、Ｉｔｅｒａｔｉｖｅ  ＳＷＩＭ法（以下、単にＳＷＩＭ法と呼ぶ）と、ＭＥＤＩ法とがある。
【００２５】
ＴＫＤ法では、点ダイポール磁場ｄ（ｒ）のフーリエ成分Ｄ（ｋ）の逆数をとり、そのマジックアングル領域近傍を０または任意の値に打ち切る。その後、計測した磁場分布δ（ｒ）のフーリエ成分Δ（ｋ）と掛け合わせ、実空間に逆フーリエ変換することで磁化率分布χ（ｒ）を算出する。
【００２６】
Ｉｔｅｒａｔｉｖｅ  ＳＷＩＭ法では、ＴＫＤ法などで算出した磁化率分布に対して、閾値処理により微細構造のみを抽出する。次に、微細構造以外の磁化率をマスク処理した微細構造磁化率分布を算出する。元の磁化率分布と微細構造磁化率分布とをｋ空間上で混合し、修正磁化率分布を算出する。修正磁化率分布に対して、閾値処理により再び微細構造を抽出する。これらの一連の処理を繰り返すことで磁化率分布を算出する。
【００２７】
ＭＥＤＩ法では、以下の式（４）に示すように、磁場と磁化率との関係式を制約条件として、正則化項と呼ばれる平滑効果を与える項を繰り返し演算により最小化することで磁化率分布χ（ｒ）を算出する。
【数４】
2016076076.tif 000006
ここで、χは、算出対象の磁化率分布ベクトル、δは、位相画像から算出した磁場分布ベクトル、Ｗは、絶対値画像を基に算出する重みベクトルを対角成分に持つ行列、Ｍは、絶対値画像の空間勾配を基にエッジ領域を０、その他の領域を１とした２値エッジマスクベクトルを対角成分に持つ行列をそれぞれ表す。また、ＦDは畳み込み積分を表す演算子、∇（ナブラ）は、３次元の空間勾配演算子、‖ａ‖nは、は任意のベクトルａのＬｎノルムを表す。また、式（４）の第２項のλは、平滑化効果と磁化率精度とのバランスをとる正則化パラメータである。
【００２８】
式（４）に示すように、ＭＥＤＩ法は正則化項として、算出途中の磁化率分布の全変動（ＴＶ：Ｔｏｔａｌ  Ｖａｒｉａｔｉｏｎ）∇χに、絶対値画像に基づくエッジマスクＭを掛け合わせたベクトルのＬ１ノルムを用いる。一般的に、任意ベクトルａのＬ１ノルム正則化付き最小二乗解は、スパース（ｓｐａｒｓｅ：疎）な解となる。したがって、ＭＥＤＩ法では、絶対値画像の変化が少ない領域が平滑化され、変化の大きいエッジ部分が保存された磁化率分布が算出される。
【００２９】
ＭＥＤＩ法を含めた正則化付き最小二乗法は、正則化パラメータの値に対応した大きさの信号値を打ち切る処理と解釈できる。そのため、正則化パラメータが大きすぎるとノイズと共に微細構造も除去され、さらに画像全体の磁化率も減少するため、算出磁化率が低下する。また、正則化パラメータが小さすぎると微細構造や画像全体の磁化率は維持するもののノイズが増加する。
【００３０】
したがって、ＭＥＤＩ法では、適切な正則化パラメータを設定する必要がある。適切な正則化パラメータを設定した場合、式（４）の第１項の残差の標準偏差は、計測した磁場分布に混入しているノイズの標準偏差に等しくなる。しかしながら、通常、磁場分布に混入したノイズの情報を事前に得ることができないため、正則化パラメータを複数回変化させて算出した残差の変化曲線から適切な正則化パラメータを算出する。従って、ＭＥＤＩ法では、適切な正則化パラメータの設定のために膨大な計算時間が必要となる。
【００３１】
本実施形態では、上記公知の各手法を踏まえ、ＭＲＩを用いて磁化率分布を算出するにあたり、簡易な処理で、アーチファクトおよびノイズを低減しつつ、高確度で高精度な磁化率分布を算出できるＭＲＩ装置およびＱＳＭ方法を提供する。
【００３２】
［ＭＲＩ装置の外観］
まず、本実施形態の磁気共鳴撮影装置（ＭＲＩ装置）について説明する。図１（ａ）〜図１（ｃ）は、本実施形態のＭＲＩ装置の外観図である。図１（ａ）は、ソレノイドコイルで静磁場を生成するトンネル型磁石を用いた水平磁場方式のＭＲＩ装置１００である。図１（ｂ）は、開放感を高めるために磁石を上下に分離したハンバーガー型（オープン型）の垂直磁場方式のＭＲＩ装置１２０である。また、図１（ｃ）は、図１（ａ）と同じトンネル型磁石を用い、磁石の奥行を短くし、かつ、斜めに傾けることによって、開放感を高めたＭＲＩ装置１３０である。
【００３３】
本実施形態では、これらの外観を有するＭＲＩ装置のいずれを用いることもできる。
【００３４】
なお、これらは一例であり、本実施形態のＭＲＩ装置はこれらの形態に限定されるものではない。本実施形態では、装置の形態やタイプを問わず、公知の各種のＭＲＩ装置を用いることができる。
【００３５】
以下、特に区別する必要がない場合は、ＭＲＩ装置１００で代表する。また、本実施形態では、例えば、ＭＲＩ装置１００の静磁場方向をｚ方向、それに垂直な２方向のうち、測定対称の被検体を載置するベッド面に平行な方向をｘ方向、他方向をｙ方向とする座標系を用いる。また、以下静磁場を単に磁場とも呼ぶ。
【００３６】
［ＭＲＩ装置の機能構成］
図２は、本実施形態のＭＲＩ装置１００の機能構成図である。本図に示すように、本実施形態のＭＲＩ装置１００は、被検体１０１が置かれる空間に静磁場を生成する、例えば、静磁場コイル、静磁場生成磁石などで構成される静磁場コイル１０２と、静磁場分布を調整するシムコイル１０４と、被検体１０１の計測領域に対し高周波磁場を送信する送信用高周波コイル１０５（以下、単に送信コイルという）と、被検体１０１から生じる核磁気共鳴信号を受信する受信用高周波コイル１０６（以下、単に受信コイルという）と、被検体１０１から生じる核磁気共鳴信号に位置情報を付加するために、ｘ方向、ｙ方向、ｚ方向それぞれに傾斜磁場を印加する傾斜磁場コイル１０３と、送信機１０７と、受信機１０８と、計算機１０９と、傾斜磁場用電源部１１２と、シム用電源部１１３と、シーケンス制御装置１１４と、を備える。
【００３７】
静磁場コイル１０２は、図１（ａ）、図１（ｂ）、図１（ｃ）にそれぞれ示した各ＭＲＩ装置１００、１２０、１３０の構造に応じて、種々の形態のものが採用される。傾斜磁場コイル１０３及びシムコイル１０４は、それぞれ傾斜磁場用電源部１１２及びシム用電源部１１３により駆動される。なお、本実施形態では、送信コイル１０５と受信コイル１０６とに別個のものを用いる場合を例にあげて説明するが、送信コイル１０５と受信コイル１０６との機能を兼用する１つのコイルで構成してもよい。送信コイル１０５が照射する高周波磁場は、送信機１０７により生成される。受信コイル１０６が検出した核磁気共鳴信号は、受信機１０８を通して計算機１０９に送られる。
【００３８】
シーケンス制御装置１１４は、傾斜磁場コイル１０３の駆動用電源である傾斜磁場用電源部１１２、シムコイル１０４の駆動用電源であるシム用電源部１１３、送信機１０７及び受信機１０８の動作を制御し、傾斜磁場、高周波磁場の印加および核磁気共鳴信号の受信のタイミングを制御する。制御のタイムチャートはパルスシーケンスと呼ばれ、計測に応じて予め設定され、後述する計算機１０９が備える記憶装置等に格納される。
【００３９】
計算機１０９は、ＭＲＩ装置１００全体の動作を制御するとともに、受信した核磁気共鳴信号に対して様々な演算処理を行う。本実施形態では、少なくとも１つ以上のエコー時間の複素画像や位相画像、磁場分布、磁化率分布などを生成する。計算機１０９は、ＣＰＵ、メモリ、記憶装置などを備える情報処理装置であり、計算機１０９にはディスプレイ１１０、外部記憶装置１１１、入力装置１１５などが接続される。
【００４０】
ディスプレイ１１０は、演算処理で得られた結果等をオペレータに表示するインタフェースである。入力装置１１５は、本実施形態で実施する計測や演算処理に必要な条件、パラメータ等をオペレータが入力するためのインタフェースである。本実施形態の入力装置１１５では、ユーザーが、計測するエコーの数や、基準のエコー時間、エコー間隔などの計測パラメータを入力できる。外部記憶装置１１１は、記憶装置とともに、計算機１０９が実行する各種の演算処理に用いられるデータ、演算処理により得られるデータ、入力された条件、パラメータ等を保持する。
【００４１】
［計算機の機能ブロック］
上述のように、本実施形態では、磁場と磁化率との関係式に基づく制約条件下で、平滑化処理を繰り返し演算することにより、簡易に高い精度で磁化率分布を算出する。次に、これを実現する本実施形態の計算機１０９の機能構成を説明する。
【００４２】
図３は、本実施形態の計算機１０９の機能ブロック図である。本実施形態の計算機１０９は、高周波磁場パルスの照射に応じて被検体から発生する核磁気共鳴信号（エコー信号）を複素信号として計測する計測制御部３１０と、計測制御部３１０が計測した複素信号から画素値が複素数である複素画像を再構成する画像再構成部３２０と、画像再構成部３２０が再構成した複素画像から位相画像を生成し、当該位相画像から磁場分布を算出する磁場分布算出部３３０と、磁場分布算出部３３０が算出した磁場分布から磁化率分布を算出する磁化率分布算出部３４０とを備える。
【００４３】
なお、計算機１０９は、ＣＰＵとメモリと記憶装置とを備える。そして、計算機１０９が実現する各種の機能は、記憶装置が保持するプログラムを、ＣＰＵがメモリにロードして実行することにより実現される。各機能の処理に用いる各種のデータ、処理中に生成される各種のデータは、記憶装置あるいは外部記憶装置１１１に格納される。
【００４４】
また、計算機１０９が実現する各種の機能のうち、少なくとも一つの機能は、ＭＲＩ装置１００とは独立した、情報処理装置であって、ＭＲＩ装置１００とデータの送受信が可能な情報処理装置により実現されていてもよい。さらに、全部または一部の機能は、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ  Ｓｐｅｃｉｆｉｃ  Ｉｎｔｅｇｒａｔｅｄ  Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ  Ｇａｔｅ  Ａｒｒａｙ）などのハードウェアによって実現してもよい。
【００４５】
［ＱＳＭ処理の流れ］
本実施形態の計算機１０９の上記各機能によるＱＳＭ法による処理（ＱＳＭ処理）の流れを説明する。図４は、本実施形態のＱＳＭ処理の流れを説明するためのフローチャートである。
【００４６】
まず、計測制御部３１０は、予め定められたパルスシーケンスにしたがってシーケンス制御装置１１４を制御し、予め定めたエコー時間のエコー信号を計測する（ステップＳ１００１）。
【００４７】
その後、画像再構成部３２０は、得られたエコー信号を、ｋ空間上に配置してフーリエ変換することにより、複素画像Ｉを再構成する（ステップＳ１００２）。
【００４８】
磁場分布算出部３３０は、複素画像Ｉから磁場分布δを算出する磁場分布算出処理を行う（ステップＳ１００３）。
【００４９】
その後、磁化率分布算出部３４０は、算出した磁場分布δから、磁化率分布χを算出する磁化率分布算出処理を行い（ステップＳ１００４）、算出した磁化率分布χをディスプレイ１１０に表示する（ステップＳ１００５）。
【００５０】
なお、磁化率分布χをディスプレイ１１０に表示する際、必要に応じて、ステップＳ１００３で算出した磁場分布δなどのほか、ＱＳＭ処理の過程によって算出される画像をディスプレイ１１０に表示してもよい。
【００５１】
以下、本実施形態のＱＳＭ処理の各処理の詳細を説明する。
【００５２】
［計測：Ｓ１００１］
計測制御部３１０は、入力装置１１５を介してユーザーが入力したパラメータに基づいて設定されるパルスシーケンスに従って、シーケンス制御装置１１４を動作させ、予め定めたエコー時間（ＴＥ）の核磁気共鳴信号（エコー信号）を得る、計測を行う。シーケンス制御装置１１４は、計測制御部３１０からの指示に従って、ＭＲＩ装置１００の各部を制御して計測を行う。本実施形態では、少なくとも１つのエコー時間のエコー信号を得る。
【００５３】
ここで、計測制御部３１０が計測に用いるパルスシーケンスの例を説明する。本実施形態では、例えば、ＧｒＥ（Ｇｒａｄｉｅｎｔ  Ｅｃｈｏ）系のパルスシーケンスを用いる。このＧｒＥ系のパルスシーケンスは、生体組織間の磁化率差によって生じる局所的な磁場変化に鋭敏である。
【００５４】
図５に、ＧｒＥ系のパルスシーケンスの一例として、ＲＳＳＧ（ＲＦ−ｓｐｏｉｌｅｄ−Ｓｔｅａｄｙ−ｓｔａｔｅ  Ａｃｑｕｉｓｉｔｉｏｎ  ｗｉｔｈ  Ｒｅｗｏｕｎｄ  Ｇｒａｄｉｅｎｔ−Ｅｃｈｏ）シーケンス５１０例を示す。本図において、ＲＦ、Ｇｓ、Ｇｐ、Ｇｒはそれぞれ、高周波磁場、スライス傾斜磁場、位相エンコード傾斜磁場、リードアウト傾斜磁場を表す。
【００５５】
ＲＳＳＧシーケンス５１０では、スライス傾斜磁場パルス５０１の印加とともに高周波磁場（ＲＦ）パルス５０２を照射し、被検体１０１内の所定のスライスの磁化を励起する。次いで磁化の位相にスライス方向および位相エンコード方向の位置情報を付加するためのスライスエンコード傾斜磁場パルス５０３および位相エンコード傾斜磁場パルス５０４を印加する。
【００５６】
画素内の核磁化の位相を分散させるディフェーズ用のリードアウト傾斜磁場パルス５０５を印加した後、リードアウト方向の位置情報を付加するためのリードアウト傾斜磁場パルス５０６を印加しながら核磁気共鳴信号（エコー）５０７を一つ計測する。そして最後に、スライスエンコード傾斜磁場パルス５０３および位相エンコード傾斜磁場パルス５０４によってディフェーズされた核磁化の位相を収束させるリフェーズ用のスライスエンコード傾斜磁場パルス５０８および位相エンコード傾斜磁場パルス５０９を印加する。
【００５７】
計測制御部３１０は、以上の手順を、スライスエンコード傾斜磁場パルス５０３、５０８（スライスエンコード数ｋｓ）および位相エンコード傾斜磁場パルス５０４、５０９（位相エンコード数ｋｐ）の強度と、ＲＦパルス５０２の位相とを変化させながら、繰り返し時間ＴＲで繰り返し実行し、１枚の画像を得るために必要なエコーを計測する。なお、このとき、ＲＦパルス５０２の位相は、例えば、１１７度ずつ増加させる。また、図５において、ハイフン以下の数字は、繰り返しの何回目であるかを示す。
【００５８】
なお、計測された各エコーはｋｒ、ｋｐ、ｋｓを座標軸とする３次元のｋ空間上に配置される。このとき、一つのエコーはｋ空間上でｋｒ軸に平行な１ラインを占める。このＲＳＳＧシーケンス５１０により得られる絶対値画像は、ＴＥを短く設定するとＴ１（縦緩和時間）強調画像、ＴＥを長く設定すると画素内の位相分散を反映したＴ２*強調画像となる。
【００５９】
なお、ここでは、ｋ空間の座標軸に平行にデータを取得するカーテシアン撮像の一つであるＲＳＳＧシーケンス５１０を例示するが、本実施形態では、任意のシーケンスを用い、任意のｋ空間領域のデータを取得することができる。
【００６０】
例えば、一回のＴＲで複数のエコーを取得し、それらを異なるｋ空間に配置するシーケンスを用いてもよい。また、計測を複数回行い、それぞれの計測で別々のｋ空間領域のデータを取得し、それらを統合し、１枚の画像を得るシーケンスであってもよい。また、ｋ空間において回転状にデータを取得するラジアルスキャンなど、ノンカーテシアン撮像を用いてもよい。
【００６１】
さらに、複数のＴＥで取得した別々の複素画像や位相画像から一枚の複素画像や位相画像を作成するシーケンスを用いてもよい。例えば、複数のＴＥで得た信号から画像を得るＤｉｘｏｎ法を実現可能なパルスシーケンスを用いる場合、水と脂肪とを分離する過程で磁場分布δが得られる。従って、この場合、後述の位相画像Ｐを生成する処理は不要となる。
【００６２】
［画像再構成：Ｓ１００２］
次に、上記ステップＳ１００２の画像再構成処理を説明する。画像再構成部３２０は、ステップＳ１００１で計測したエコー時間のエコー信号をｋ空間上に配置し、フーリエ変換する。これにより、画像再構成部３２０は、各画素値が複素数の複素画像Ｉを算出する。
【００６３】
［磁場分布算出処理：Ｓ１００３］
次に、上記ステップＳ１００３の、磁場分布算出部３３０による磁場分布算出処理を説明する。本実施形態の磁場分布算出部３３０は、画像再構成部３２０が再構成した複素画像Ｉから、磁場分布δを算出する。算出される磁場分布δは、生体組織間の磁化率差によって生じる局所的な磁場変化をとらえたものである。
【００６４】
磁場分布算出部３３０は、まず、複素画像Ｉから位相画像Ｐを算出し、算出した位相画像Ｐに、背景位相を除去するための各種の位相画像処理を行う。そして、位相画像処理後の位相画像Ｐから、磁場分布δを得る。
【００６５】
位相画像処理では、磁化率以外に由来する位相成分を除き、磁化率に由来する位相成分を抽出する。具体的には、算出した位相画像Ｐにおいて、−πからπの範囲を超えて折り返された位相を除去する位相アンラップ処理と、生体形状に起因して生じる大域的な位相変化を除去する大域的位相変化除去処理と、を行う。
【００６６】
これらを実現するため、本実施形態の磁場分布算出部３３０は、図３に示すように、複素画像Ｉから位相画像Ｐを算出する位相画像算出部３３１と、位相アンラップ処理を行う位相アンラップ処理部３３２と、大域的位相変化除去処理を行う大域的位相変化除去部３３３と、位相処理後の位相画像Ｐを磁場分布δに変換する磁場分布変換部３３４と、を備える。
【００６７】
位相画像算出部３３１は、複素画像Ｉの各画素の画素値の偏角を求め、それを各々画素値とすることにより位相画像Ｐを算出する。すなわち、位相画像Ｐは、複素画像Ｉの各画素の複素数の位相成分を画素値とする画像である。ここで、画素ｒにおける位相画像Ｐの画素値Ｐ（ｒ）は、複素画像Ｉの画素値Ｉ（ｒ）を用いて、式（５）により算出される。
Ｐ（ｒ）＝ａｒｇ｛Ｉ（ｒ）｝  ・・・（５）
【００６８】
位相アンラップ処理部３３２は、位相画像Ｐに対し、−πからπの範囲を超えて折り返された位相を除去するための位相アンラップ処理を実行する。位相画像Ｐにおける一部の領域では、−πからπの範囲を超えた位相値が−πからπの範囲内に折り返される。そのため、撮像部位（例えば、頭部など）の全領域において正確な位相を求めるためには、これらの折り返しを補正する必要がある。本実施形態では、例えば領域拡大法などを用い、−πからπの範囲に折り返されている位相値を補正する。
【００６９】
大域的位相変化除去部３３３は、位相アンラップ処理後の位相画像Ｐにおいて、生体形状に起因して生じる大域的な位相変化を除去する。
【００７０】
位相画像Ｐの各画素値は、撮像部位の形状等に依存して生じる静磁場不均一に起因した大域的な位相変化と、組織間の磁化率変化に起因する局所的な位相変化との和である。大域的な位相変化と局所的な位相変化とは、それぞれ、位相画像Ｐのｋ空間上における低周波成分と高周波成分に対応する。本処理によって、大域的な位相変化のみを除去し、組織間の磁化率変化に起因する局所的な位相変化を抽出する。
【００７１】
本実施形態の大域的位相変化除去部３３３は、大域的位相変化除去処理として、まず、撮像した三次元画像（元画像）に対して二次元画像ごとにそれぞれローパスフィルタ処理を施し、低解像度画像を算出する。そして、元画像を低解像度画像で複素除算する。これにより、低解像度画像に含まれる大域的な位相変化を元画像から除去する。
【００７２】
なお、大域的な位相変化を除去する方法には様々な方法がある。例えば、三次元画像を低次多項式でフィッティングすることにより大域的な位相変化を抽出し、それを元画像から減算する方法などがある。大域的位相変化除去では、こうした他の方法を用いてもよい。
【００７３】
なお、位相画像処理の、位相アンラップ処理と、大域的位相変化除去処理は、あくまで一例であり、これらに限定されない。また、これら２つの処理のうちのどちらかの処理を省略することもできる。また、各処理の処理順は問わない。また、これら２つの処理は、行わなくてもよい。
【００７４】
磁場分布変換部３３４は、式（１）を用いて、位相画像処理後の位相画像Ｐを、磁場分布δに変換する。
【００７５】
［磁場分布算出処理の流れ］
本実施形態の磁場分布算出部３３０の各機能による磁場分布算出処理の流れの一例を説明する。図６は、本実施形態の磁場分布算出処理の処理フローである。
【００７６】
位相画像算出部３３１は、複素画像Ｉから位相画像Ｐを算出する（ステップＳ１１０１）。
【００７７】
そして、磁場分布算出部３３０は、位相画像Ｐに、各種の位相画像処理を行い、位相処理後の位相画像を得る。本実施形態では、例えば、位相アンラップ処理部３３２が、算出した位相画像Ｐに対し、位相アンラップ処理を施し（ステップＳ１１０２）、大域的位相変化除去部３３３が、位相アンラップ処理後の位相画像Ｐに対し、大域的位相変化除去処理を行う（ステップＳ１１０３）。上述のように、これらの位相画像処理の処理順は問わない。また、位相画像処理は、いずれかのみを行ってもよいし、全て行わなくてもよい。さらに、これらに限定されない。
【００７８】
最後に、磁場分布変換部３３４は、位相処理後の位相画像Ｐを磁場分布に変換し（ステップＳ１１０４）、処理を終了する。
【００７９】
［磁化率分布算出処理：Ｓ１００４］
次に、上記ステップＳ１００４の、磁化率分布算出部３４０による、磁化率分布算出処理を説明する。本実施形態の磁化率分布算出部３４０は、磁場分布算出部３３０で算出した磁場分布δを用いて磁化率分布χを算出する。
【００８０】
磁化率分布算出部３４０は、予め定めた初期磁化率分布χ0と、磁場分布δとにより定義した誤差関数を繰り返し演算により最小化する最小化処理により磁化率分布χを算出する。このとき、最小化処理中に、誤差関数に基づいて算出される更新磁化率分布χUPDを、当該更新磁化率分布が算出される毎に平滑化する。なお、更新磁化率分布χUPDは、初期磁化率分布χ0を誤差関数に基づいて更新した磁化率分布である。
【００８１】
これを実現するため、本実施形態の磁化率分布算出部３４０は、誤差関数を設定する誤差関数設定部３４１と、設定された誤差関数を用いて、磁場分布δと、初期磁化率分布χ0とから、更新磁化率分布χUPDを算出する更新磁化率分布算出部３４２と、更新磁化率分布χUPDを平滑化し、合成磁化率分布χADDを生成する平滑化処理部３４３と、繰り返し演算の終了判定を行うとともに、初期磁化率分布χ0を算出した合成磁化率分布χADDに更新する判定部３４６と、を備える。
【００８２】
また、平滑化処理部３４３は、更新磁化率分布χUPDに、予め定めた平滑化フィルタを適用して、フィルタ処理後の磁化率分布（処理後磁化率分布）χSMTを生成するフィルタ処理部３４４と、更新磁化率分布χUPDと、処理後磁化率分布χSMTとを合成し、合成磁化率分布χADDを生成する合成部３４５と、を備える。
【００８３】
［誤差関数設定部］
式（２）に示すように、磁場分布δは、磁化率分布χを空間的に畳み込み積分したものである。位相画像Ｐ内の全ての画素を対象とするため、式（２）をベクトルと行列とにより表現すると、以下の式（６）のように表される。
δ＝ＦDχ  ・・・（６）
ここで、δは、全画素数Ｎの大きさを持つ磁場分布の列ベクトル、χは、算出する磁化率分布の列ベクトルである、また、ＦDは、Ｎ×Ｎの大きさを持ち、磁化率分布χに対する畳み込み演算に相当する行列である。
【００８４】
本実施形態の更新磁化率分布算出部３４２は、式（６）から最小二乗法によって磁化率分布χを求める。このとき、誤差関数を導入し、それを最小化することにより、解を得る。
【００８５】
本実施形態の誤差関数設定部３４１は、更新磁化率分布算出部３４２が更新磁化率分布算出の際に用いる誤差関数を設定する。設定される誤差関数ｅ（χ0）は、例えば、以下の式（７）で表されるものとする。
【数７】
2016076076.tif 000007
ここで、Ｗは、複素画像Ｉの位相値成分を基に算出される重みベクトルを対角成分に持つ行列（以下、重みＷと呼ぶ）、ＦDは、畳み込み積分を表す演算子、‖ａ‖22は、は任意のベクトルａのＬ２ノルム（ユークリッドノルム）の二乗をそれぞれ表す。
【００８６】
すなわち、誤差関数は、磁場分布δと初期磁化率分布χ0に畳み込み演算を行った結果とにより定まる関数に、重みＷを付与したものである。この重みＷは、上述のように複素画像の位相値成分（位相画像Ｐ）に基づいて算出される。
【００８７】
具体的には、誤差関数設定部３４１は、位相画像Ｐから、各画素の位相ばらつきを算出し、算出した位相ばらつきの大きさに基づいて重みＷを算出する。例えば、位相ばらつきが大きいほど、重みＷが小さくなるよう算出する。
【００８８】
位相ばらつきには、位相ばらつきを算出する対象画素の周辺領域の、複数の画素の画素値（位相値）を用いて算出した、例えば、標準偏差等の統計指標を用いる。
【００８９】
なお、用いる位相画像Ｐは、複素画像Ｉから得たものであってもよいし、複素画像Ｉから得た位相画像Ｐに、上述の各種の位相画像処理を施したものであってもよい。
【００９０】
［更新磁化率算出部］
更新磁化率分布算出部３４２は、上述のように、誤差関数ｅ（χ0）を用いて、初期磁化率分布χ0から更新磁化率分布χUPDを算出する。ここでは、更新磁化率分布算出部３４２は、誤差関数ｅ（χ0）の勾配∇ｅ（χ0）を算出し、算出した勾配∇ｅ（χ0）を用いて、以下の式（８）により、更新磁化率分布χUPDを算出する。
χUPD＝χ0＋∇ｅ（χ0）  ・・・（８）
【００９１】
なお、本実施形態では、最初に設定する初期磁化率分布χ0は、任意のものを設定できる。例えば、全ての値を０としたものを初期磁化率分布χ0としてもよい。
【００９２】
［平滑化処理部；フィルタ処理部］
平滑化処理部３４３は、更新磁化率分布χUPDに、平滑化フィルタを適用して平滑化し、処理後磁化率分布χSMTを得る。この中で、フィルタ処理部３４４は、平滑化フィルタとして、エッジ保存平滑化フィルタを用い、更新磁化率分布χUPDにおける、各画素値の局所的な磁化率変動に合わせて適応的に線形平滑化処理を行う。
【００９３】
具体的には、更新磁化率分布χUPDにおける任意の画素位置をｒとしたとき、その周辺のｎ×ｎ（ｎは１以上の整数）の局所領域における平均値をμ（ｒ）、分散をσ（ｒ）2とすると、位置ｒにおけるエッジ保存平滑化フィルタ適用後の磁化率分布である処理後磁化率分布χSMTの画素値χSMT（ｒ）は、以下の式（９）で表される。
【数９】
2016076076.tif 000008
ここでν2は、磁化率分布χに混入したノイズ分散に相当するパラメータである。本実施形態では、このパラメータとして、更新磁化率分布χUPDの各画素における局所分散σ（ｒ）2の平均値を採用する。
【００９４】
本実施形態で適用するエッジ保存平滑化フィルタによれば、式（９）に示すように、局所分散が小さい領域では平滑化が強く働き、局所平均値に近い値となるため、磁化率の定量値を維持できる。また、局所分散が大きい領域では平滑化が弱く働き、微細構造を含むエッジ部分を保存する。したがって、更新磁化率分布χUPDにこのエッジ保存平滑化フィルタを適用することにより、磁化率分布のノイズを低減し、定量値を維持しながら、微細静脈や磁化率変化の大きい組織間のエッジを自動的に抽出し保存することができる。
【００９５】
なお、本実施形態では、上述のエッジ保存平滑化フィルタを用いたが、平滑化処理に用いるフィルタは、これに限定されない。例えばメディアンフィルタやバイラテラルフィルタなど公知のエッジ保存平滑化フィルタを使用してもよい。
【００９６】
［平滑化処理部；合成部］
合成部３４５は、更新磁化率分布χUPDと処理後磁化率分布χSMTとを合成して合成磁化率分布χADDを算出する。合成は、ノイズの増加を防ぎ、かつ、磁化率分布の微細構造を維持するよう行う。例えば、ノイズの影響が大きい、ｋ空間のマジックアングル近傍にエッジ保存平滑化された処理後磁化率分布χSMTのフーリエ成分が挿入されるよう、合成する。
【００９７】
合成部３４５は、更新磁化率分布χUPDおよび処理後磁化率分布χSMTをｋ空間データに変換し、変換後の更新磁化率分布χUPDおよび処理後磁化率分布χSMTのｋ空間データそれぞれに重みをつけて加算し、合成磁化率分布χADDを得る。
【００９８】
具体的には、合成部３４５は、まず、更新磁化率分布χUPDおよび処理後磁化率分布χSMTそれぞれを、フーリエ変換によりｋ空間データに変換する。そして、変換後の各ｋ空間データ（フーリエ成分）に、上記目的を達成するよう重み付けを行い、重み付けを行った各ｋ空間データを加算する。そして、加算後、逆フーリエ変換により合成磁化率分布χADDを得る。
【００９９】
ｋ空間における重み付けは、重み画像Ｇを作成し、重み画像Ｇをｋ空間データに掛け合わせることにより行う。このとき、更新磁化率分布χUPDから得たｋ空間データに乗算する重み画像ＧUPDでは、例えば、マジックアングル領域近傍のデータに付与する重みは、他の領域のデータに付与する重みより小さく設定される。また、処理後磁化率分布χSMTから得たｋ空間データに乗算する重み画像ＧSMTでは、逆に、マジックアングル領域近傍のデータに付与する重みは、他の領域のデータに付与する重みより大きく設定される。
【０１００】
これは、更新磁化率分布χUPDは、比較的ノイズが多いため、マジックアングル領域近傍のデータの寄与を低減させ、ノイズの増加を防ぐためである。また、処理後磁化率分布χSMTは、ノイズが低減されているため、マジックアングル領域近傍のデータの寄与を増大させ、微細構造を維持するためである。
【０１０１】
これらの条件を満たす重み画像ＧUPD（ｋ）、ＧSMT（ｋ）として、点ダイポール磁場のフーリエ成分Ｄ（ｋ）とマジックアングル領域の大きさを調整する調整パラメータｂとを用いて、例えば、以下の式（１０）より定義する重み画像を用いることができる。なお、ｋは、ｋ空間のデータ位置を示す。
【数１０】
2016076076.tif 000009
上記式（１０）では、更新磁化率分布χUPDにおいて、フーリエ成分の値が所定の値（ｂ）より小さい領域を、マジックアングル領域とする。
【０１０２】
なお、調整パラメータｂが０の場合は、ＧUPD（ｋ）が１となり、ＧSMT（ｋ）が０となる。このため、エッジ保存平滑化処理を実施せずに式（６）の誤差関数を最小化する単純な最小二乗法となる。
【０１０３】
［判定部］
判定部３４６は、繰り返し演算により誤差関数ｅ（χ0）を最小化する最小化処理において、繰り返し演算の終了判定を行う。処理を継続する、すなわち、繰り返すと判定した場合は、その時点で得られた最新の合成磁化率分布χADDを、初期磁化率分布χ0として設定することにより、初期磁化率分布χ0を更新する。なお、終了と判定した際は、その時点の最新の合成磁化率分布χADDを磁化率分布として出力する。
【０１０４】
本実施形態では、判定部３４６は、更新磁化率分布χUPDが算出されると、その更新磁化率分布χUPDとその時点の初期磁化率分布χ0とを用いて、終了か否かを判定する。判定部３４６は、例えば、以下の式（１１）に示す評価関数ｆ（χUPD）を用い、その値が予め定めた閾値εより小さい場合、終了と判定する。
【数１１】
2016076076.tif 000010
【０１０５】
［磁化率分布算出処理の流れ］
次に、本実施形態の磁化率分布算出部３４０による、磁化率分布算出処理の流れを説明する。図７は、本実施形態の磁化率分布算出処理の処理フローである。また、図８は、磁化率分布算出処理を説明するための説明図である。
【０１０６】
本実施形態では、予め定めた初期磁化率分布χ0と磁場分布δとにより定義した誤差関数を繰り返し演算により最小化する最小化処理を行い、磁化率分布χを算出する。このとき、最小化処理中に、誤差関数から算出される更新磁化率分布に対して平滑化処理を行い、合成磁化率分布を得、この合成磁化率分布で初期磁化率分布を更新し、処理を繰り返す。具体的には、以下の処理を行う。
【０１０７】
磁化率分布算出部３４０は、初期磁化率分布χ0４１０を設定する（ステップＳ１２０１）。そして、誤差関数設定部３４１は、初期磁化率分布χ0４１０と磁場分布δ４００とを用いて誤差関数ｅ（χ0）を設定する（ステップＳ１２０２）。
【０１０８】
そして、更新磁化率分布算出部３４２は、誤差関数の勾配∇ｅ（χ0）を算出し（ステップＳ１２０３）、上記式（８）に従って、初期磁化率分布χ0を更新し、更新磁化率分布χUPD４２０を算出する（ステップＳ１２０４）。
【０１０９】
判定部３４６は、上記式（１１）の評価関数ｆ（χUPD）を用いて、終了判定を行う（ステップＳ１２０５）。
【０１１０】
ステップＳ１２０５において、継続するものと判定された場合、フィルタ処理部３４４は、更新磁化率分布χUPD４２０に、予め定めた平滑化フィルタを適用し（Ｆｉｌｔｅｒｉｎｇ）、処理後磁化率分布χSMT４３０を算出する（ステップＳ１２０６）。
【０１１１】
合成部３４５は、更新磁化率分布χUPD４２０と処理後磁化率分布χSMT４３０とを合成し、合成磁化率分布χADD４４０を得る（ステップＳ１２０７）。合成は、更新磁化率分布χUPD４２０および処理後磁化率分布χSMT４３０をそれぞれフーリエ変換（ＦＴ）して得たｋ空間データ４２１、４３１に、上記式（１０）の重み画像４２２、４３２を乗算して加算し、その重み付け加算の結果を逆フーリエ変換（ＩＦＴ）することにより行う。
【０１１２】
判定部３４６は、合成磁化率分布χADD４４０を初期磁化率分布χ0４１０として設定し（ステップＳ１２０８）、ステップＳ１２０２に戻り、処理を繰り返す。
【０１１３】
ここで、ステップＳ１２０５において、終了と判定された場合は、判定部３４６は、その時点で得られている最新の合成磁化率分布χADDを、磁化率分布算出部３４０が出力する磁化率分布χとし（ステップＳ１２０９）、処理を終了する。
【０１１４】
なお、ステップＳ１２０５において、合成磁化率分布χADDが算出されていない場合、すなわち、１回目の演算で、終了と判定された場合は、その時点の更新磁化率分布χUPDから合成磁化率分布χADDを上記手法で算出し、算出した合成磁化率分布χADDを、出力する磁化率分布χとする。
【０１１５】
＜磁化率分布算出処理の流れの変形例＞
なお、判定部３４６は、合成磁化率分布χADDを用いて終了か否かを判定してもよい。すなわち、合成磁化率分布χADDが算出されると、当該合成磁化率分布χADDとその時点の初期磁化率分布χ0とを用いて、終了するか否かを判定してもよい。
【０１１６】
この場合、終了を判定する評価関数として、以下のｆ（χADD）を用いる。
【数１２】
2016076076.tif 000011
【０１１７】
この場合の、磁化率分布算出部３４０の各部による磁化率分布算出処理の流れを、図９に示す。また、図１０は、本処理を説明するための説明図である。
【０１１８】
磁化率分布算出部３４０は、初期磁化率分布χ0４１０を設定する（ステップＳ１３０１）。そして、誤差関数設定部３４１は、初期磁化率分布χ0４１０と磁場分布δ４００とを用いて誤差関数ｅ（χ0）を設定する（ステップＳ１３０２）。
【０１１９】
そして、更新磁化率分布算出部３４２は、誤差関数の勾配∇ｅ（χ0）を算出し（ステップＳ１３０３）、上記式（８）に従って、更新磁化率分布χUPD４２０を算出する（ステップＳ１３０４）。
【０１２０】
フィルタ処理部３４４は、更新磁化率分布χUPD４２０に、予め定めた平滑化フィルタを適用し（Ｆｉｌｔｅｒｉｎｇ）、処理後磁化率分布χSMT４３０を算出する（ステップＳ１３０５）。
【０１２１】
合成部３４５は、更新磁化率分布χUPD４２０と処理後磁化率分布χSMT４３０とを合成し、合成磁化率分布χADD４４０を得る（ステップＳ１３０６）。合成は、更新磁化率分布χUPD４２０および処理後磁化率分布χSMT４３０をそれぞれフーリエ変換（ＦＴ）して得たｋ空間データ４２１、４３１に、上記式（１０）の重み画像４２２、４３２を乗算して加算し、その重み付け加算の結果を逆フーリエ変換（ＩＦＴ）することにより行う。
【０１２２】
次に、判定部３４６は、上記式（１２）の評価関数ｆ（χADD）を用いて、終了判定を行う（ステップＳ１３０７）。
【０１２３】
ステップＳ１３０７において、処理継続と判定された場合は、磁化率分布算出部３４０は、合成磁化率分布χADD４４０を初期磁化率分布χ0４１０として設定し（ステップＳ１３０８）、ステップＳ１３０２に戻り、処理を繰り返す。
【０１２４】
ここで、ステップＳ１３０７において、終了と判定された場合は、その時点で得られている最新の合成磁化率分布χADDを、磁化率分布算出部３４０が出力する磁化率分布χとし（ステップＳ１３０９）、処理を終了する。
【０１２５】
合成磁化率分布χADDで判定を行うと、更新磁化率分布χUPDで判定を行う場合に比べ、フィルタ処理の効果の大きい画像を得ることができる。一方更新磁化率分布χUPDで判定を行うと、合成磁化率分布χADDで判定を行う場合に比べ、実際の磁化率分布に近いデータで判定がなされるため、より適切な範囲で判定が行われる。
【０１２６】
以上説明したように、本実施形態のＭＲＩ装置は、高周波磁場パルスの照射に応じて被検体１０１から発生する核磁気共鳴信号を複素信号として計測する計測制御部３１０と、前記複素信号から画素値が複素数である複素画像Ｉを再構成する画像再構成部３２０と、前記複素画像Ｉから位相画像Ｐを生成し、当該位相画像Ｐから磁場分布δを算出する磁場分布算出部３３０と、前記磁場分布δから磁化率分布χを算出する磁化率分布算出部３４０と、を備え、前記磁化率分布算出部３４０は、予め定めた初期磁化率分布χ0と前記磁場分布δとを用いて定義した誤差関数ｅ（χ0）を繰り返し演算により最小化する最小化処理を行って前記磁化率分布χを算出するとともに、当該最小化処理中に前記誤差関数ｅ（χ0）に基づいて算出される更新磁化率分布χUPDを、当該更新磁化率分布χUPDが算出される毎に平滑化する。
【０１２７】
前記磁化率分布算出部３４０は、前記誤差関数ｅ（χ0）を用いて前記初期磁化率分布χ0から前記更新磁化率分布χUPDを算出する更新磁化率分布３４２と、前記更新磁化率分布χUPDを平滑化して合成磁化率分布χADDを生成する平滑化処理部３４３と、前記繰返し演算の終了判定を行うとともに、前記合成磁化率分布χADDを前記初期磁化率分布χ0に設定する判定部３４６と、を備え、前記平滑化処理部３４３は、前記更新磁化率分布χUPDに予め定めた平滑化フィルタを適用し、処理後磁化率分布χSMTを得るフィルタ処理部３４４と、前記更新磁化率分布χUPDと前記処理後磁化率分布χSMTとを合成し、合成磁化率分布χADDを生成する合成部３４５と、を備え、前記判定部３４６は、終了と判定した際の前記合成磁化率分布χADDを、磁化率分布χとしてもよい。
【０１２８】
さらに、前記誤差関数ｅ（χ0）は、前記磁場分布δと前記初期磁化率分布χ0に畳み込み演算を行った結果とにより定まる関数に重みＷを付与して設定され、前記重みＷは、前記複素画像Ｉから得られる位相画像Ｐに基づいて算出されてもよい。
【０１２９】
このように、本実施形態では、磁場と磁化率の関係式に基づく制約条件下で平滑化処理を繰り返し演算することにより高精度な磁化率分布を算出する。
【０１３０】
すなわち、本実施形態では、誤差関数を最小化する最小化処理の過程で、得られる更新磁化率分布χUPDを、平滑化し、処理を繰り返す。平滑化においては、平滑化フィルタにより平滑化した処理後磁化率分布χSMTと、誤差関数の勾配を用いて算出した更新磁化率分布χUPDとを、ノイズを低減し、かつ、微細構造を維持できるよう、ｋ空間上で重み付けをして合成する。
【０１３１】
このように、本実施形態によれば、平滑化フィルタによるフィルタリングにより、ノイズ部分のデータを平均化するため、算出する磁化率分布χにおけるノイズを抑えることができる。さらに、合成処理により、同磁化率分布χにおいて微細構造を維持できる。
【０１３２】
本実施形態によれば、従来手法のように、正則化項を用いないため、簡易な処理で、ノイズを抑え、かつ、微細構造が維持された磁化率分布χを得ることができる。すなわち、本実施形態によれば、適切な正則化パラメータを設定するための膨大な計算も不要である。
【０１３３】
また、本実施形態によれば、不良条件問題を解くにあたり、マジックアングル領域近傍の値を打ち切ることがない。従って、打ち切ることにより発生するアーチファクトも抑えられ、算出される磁化率分布の精度も向上する。
【０１３４】
さらに、更新磁化率分布χUPDを、ノイズを抑え、かつ、微細構造が維持されたものとするために用いる処理後磁化率分布χSMT、および、繰り返し演算の際、初期磁化率分布χ0を更新する合成磁化率分布χADDは、ともに、更新磁化率分布χUPDから算出されるものである。すなわち、位相情報に依存するデータである。
【０１３５】
本実施形態では、このような位相情報に依存するデータを用いて最小化処理の繰り返し演算中に算出される磁化率分布を更新するため、計測された位相情報と算出される磁化率分布との不整合が生じにくいという利点がある。
【０１３６】
さらに、本実施形態によれば、誤差関数内の重みＷを、位相画像から算出する。位相画像は、絶対値画像に比べ、位相ばらつきの反映精度が高いため、この位相ばらつきを精度よく反映した重みＷを得ることができる。そして、それを用いた誤差関数により磁化率分布χを算出するため、より高い精度で磁化率分布を算出することができる。
【０１３７】
以上説明したように、本実施形態によれば、簡易な処理で、磁場と磁化率の関係性を維持しながら微細構造などのエッジ情報が保存された平滑化処理を繰り返すことにより、アーチファクトおよびノイズを低減し、かつ、算出する磁化率分布の確度および精度を向上させることができる。
【０１３８】
＜合成時重み画像の変形例＞
なお、合成部３４５が更新磁化率分布χUPDと処理後磁化率分布χSMTとを合成する際に用いる重み画像Ｇは、上記の重み画像に限定されない。例えば、ＧUPD（ｋ）およびＧSMT（ｋ）を、それぞれ、一定の値としてもよい。すなわち、各ｋ空間データ内の全データに、一様な重みを設定してもよい。
【０１３９】
一様な重みを設定することにより、最終的に得られる画像のＳＮＲを制御しやすくなる。
【０１４０】
＜誤差関数内の重みの変形例＞
また、誤差関数設定部３４１が設定する誤差関数ｅ（χ0）内の重みＷは、上記実施形態では、複素画像Ｉの位相値成分から算出するものとしているが、これに限定されない。例えば、複素画像Ｉの絶対値成分から算出してもよい。絶対値成分から算出することにより、簡単に少ない処理で重みを算出することができる。
【符号の説明】
【０１４１】
１００：ＭＲＩ装置、１０１：被検体、１０２：静磁場コイル、１０３：傾斜磁場コイル、１０４：シムコイル、１０５：送信コイル、１０６：受信コイル、１０７：送信機、１０８：受信機、１０９：計算機、１１０：ディスプレイ、１１１：外部記憶装置、１１２：傾斜磁場用電源部、１１３：シム用電源部、１１４：シーケンス制御装置、１１５：入力装置、１２０：ＭＲＩ装置、１３０：ＭＲＩ装置、３１０：計測制御部、３２０：画像再構成部、３３０：磁場分布算出部、３３１：位相画像算出部、３３２：位相アンラップ処理部、３３３：大域的位相変化除去部、３３４：磁場分布変換部、３４０：磁化率分布算出部、３４１：誤差関数設定部、３４２：更新磁化率分布算出部、３４３：平滑化処理部、３４４：フィルタ処理部、３４５：合成部、３４６：判定部、４００：磁場分布、４１０：初期磁化率分布、４２０：更新磁化率分布、４２１：更新磁化率分布のｋ空間データ、４２２：更新磁化率分布に乗算する重み画像、４３０：処理後磁化率分布、４３１：処理後磁化率分布のｋ空間データ、４３２：処理後磁化率分布に乗算する重み画像、４４０：合成磁化率分布、５０１：スライス傾斜磁場パルス、５０２：ＲＦパルス、５０３：スライスエンコード傾斜磁場パルス、５０４：位相エンコード傾斜磁場パルス、５０５：リードアウト傾斜磁場パルス、５０６：リードアウト傾斜磁場パルス、５０７：エコー、５０８：スライスエンコード傾斜磁場パルス、５０９：位相エンコード傾斜磁場パルス、５１０：ＲＳＳＧシーケンス
【図１】
2016076076.tif 000012
【図２】
2016076076.tif 000013
【図３】
2016076076.tif 000014
【図４】
2016076076.tif 000015
【図５】
2016076076.tif 000016
【図６】
2016076076.tif 000017
【図７】
2016076076.tif 000018
【図８】
2016076076.tif 000019
【図９】
2016076076.tif 000020
【図１０】
2016076076.tif 000021
【手続補正書】
【提出日】平成29年4月28日(2017.4.28)
【手続補正２】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
高周波磁場パルスの照射に応じて被検体から発生する核磁気共鳴信号を複素信号として計測する計測制御部と、
前記複素信号から画素値が複素数である複素画像を再構成する画像再構成部と、
前記複素画像から位相画像を生成し、当該位相画像から磁場分布を算出する磁場分布算出部と、
前記磁場分布から磁化率分布を算出する磁化率分布算出部と、を備え、
前記磁化率分布算出部は、予め定めた初期磁化率分布と前記磁場分布とを用いて定義した誤差関数を繰り返し演算により最小化する最小化処理を行って前記磁化率分布を算出するとともに、当該最小化処理中に前記誤差関数に基づいて算出される更新磁化率分布を、
当該更新磁化率分布が算出される毎に平滑化すること
を特徴とする磁気共鳴イメージング装置。
【請求項２】
請求項１記載の磁気共鳴イメージング装置であって、
前記磁化率分布算出部は、
前記誤差関数を用いて前記初期磁化率分布から前記更新磁化率分布を算出する更新磁化率分布算出部と、
前記更新磁化率分布に平滑化フィルタを適用して平滑化し、合成磁化率分布を生成する平滑化処理部と、
前記繰返し演算の終了判定を行うとともに、前記初期磁化率分布を前記合成磁化率分布に更新する判定部と、を備え、
前記判定部は、終了と判定した際の前記合成磁化率分布を、前記磁化率分布とすることを特徴とする磁気共鳴イメージング装置。
【請求項３】
請求項２記載の磁気共鳴イメージング装置であって、
前記平滑化処理部は、
前記更新磁化率分布に予め定めた前記平滑フィルタを適用し、処理後磁化率分布を生成するフィルタ処理部と、
前記更新磁化率分布と前記処理後磁化率分布とを合成し、前記合成磁化率分布を生成する合成部を備えることを特徴とする磁気共鳴イメージング装置。
【請求項４】
請求項１乃至３のいずれかに記載の磁気共鳴イメージング装置であって、
前記誤差関数は、前記磁場分布と前記初期磁化率分布に畳み込み演算を行った結果とにより定まる関数に重みを付与して設定され、
前記重みは、前記複素画像から得られる位相画像に基づいて算出されること
を特徴とする磁気共鳴イメージング装置。
【請求項５】
請求項１乃至３のいずれかに記載の磁気共鳴イメージング装置であって、
前記磁化率分布算出部は、前記誤差関数の勾配を算出し、算出した勾配を用いて、
前記更新磁化率分布を算出すること
を特徴とする磁気共鳴イメージング装置。
【請求項６】
請求項３に記載の磁気共鳴イメージング装置であって、
前記合成部は、前記更新磁化率分布および前記処理後磁化率分布をそれぞれｋ空間データに変換し、前記変換後のｋ空間データそれぞれに重みを付けて加算することにより、前記合成磁化率分布を生成すること
を特徴とする磁気共鳴イメージング装置。
【請求項７】
請求項６記載の磁気共鳴イメージング装置であって、
前記更新磁化率分布のｋ空間データにおいて、マジックアングル領域近傍のデータに付与する重みは、他の領域のデータに付与する重みより小さく設定され、
前記処理後磁化率分布のｋ空間データにおいて、マジックアングル領域近傍のデータに付与する重みは、他の領域のデータに付与する重みより大きく設定されること
を特徴とする磁気共鳴イメージング装置。
【請求項８】
請求項６記載の磁気共鳴イメージング装置であって、
前記重みは、１つの前記ｋ空間データ内の全データに一様に設定されること
を特徴とする磁気共鳴イメージング装置。
【請求項９】
請求項２記載の磁気共鳴イメージング装置であって、
前記平滑化フィルタは、エッジ保存平滑化フィルタであること
を特徴とする磁気共鳴イメージング装置。
【請求項１０】
請求項１記載の磁気共鳴イメージング装置であって、
前記初期磁化率分布は、全ての値が０であること
を特徴とする磁気共鳴イメージング装置。
【請求項１１】
請求項２記載の磁気共鳴イメージング装置であって、
前記判定部は、前記更新磁化率分布が算出されると、当該更新磁化率分布と前記初期磁化率分布とを用いて前記終了判定を行うこと
を特徴とする磁気共鳴イメージング装置。
【請求項１２】
請求項２記載の磁気共鳴イメージング装置であって、
前記判定部は、前記合成磁化率分布が算出されると、当該合成磁化率分布と前記初期磁化率分布とを用いて前記終了判定を行うこと
を特徴とする磁気共鳴イメージング装置。
【請求項１３】
高周波磁場パルスの照射に応じて被検体から発生する核磁気共鳴信号を収集し、再構成することにより得た複素画像である磁気共鳴画像から位相画像を得、
前記位相画像から磁場分布を算出し、
予め定めた初期磁化率分布と前記磁場分布とにより定義した誤差関数を、繰り返し演算により最小化する最小化処理を行い、磁化率分布を算出し、
前記最小化処理中に、
前記誤差関数から算出される、前記初期磁化率分布を更新した更新磁化率分布に対して平滑化を行うこと
を特徴とする磁気共鳴イメージング装置における定量的磁化率マッピング方法。
【請求項１４】
磁気共鳴イメージング装置が収集した核磁気共鳴信号或いは当該核磁気共鳴信号を再構成した複素画像である磁気共鳴画像を用いて磁化率分布を算出する計算機であって、
前記複素画像から位相画像を生成し、当該位相画像から磁場分布を算出する磁場分布算出部と、
予め定めた初期磁化率分布と前記磁場分布とを用いて定義した誤差関数を繰り返し演算により最小化する最小化処理を行って前記磁化率分布を算出する磁化率分布算出部と、
前記最小化処理中に前記誤差関数に基づいて算出される更新磁化率分布を、当該更新磁化率分布が算出される毎に平滑化する平滑処理部と、
を備えることを特徴とする計算機。
【請求項１５】
磁気共鳴イメージング装置が収集した核磁気共鳴信号或いは当該核磁気共鳴信号を再構成した複素画像である磁気共鳴画像を用いて磁化率分布を算出する磁化率分布計算方法であって、
前記複素画像から位相画像を生成し、当該位相画像から磁場分布を算出するステップと、
予め定めた初期磁化率分布と前記磁場分布とを用いて定義した誤差関数を繰り返し演算により最小化する最小化処理を行って前記磁化率分布を算出するステップと、
前記最小化処理中に前記誤差関数に基づいて算出される更新磁化率分布を、当該更新磁化率分布が算出される毎に平滑化するステップと、
を含む磁化率分布計算方法。
【請求項１６】
コンピュータに、
磁気共鳴イメージング装置が収集した核磁気共鳴信号或いは当該核磁気共鳴信号を再構成した複素画像である磁気共鳴画像を用いて磁化率分布を算出するステップを実行させる磁化率分布計算プログラムであって、
前記複素画像から位相画像を生成し、当該位相画像から磁場分布を算出するステップと、
予め定めた初期磁化率分布と前記磁場分布とを用いて定義した誤差関数を繰り返し演算により最小化する最小化処理を行って前記磁化率分布を算出するステップと、
前記最小化処理中に前記誤差関数に基づいて算出される更新磁化率分布を、当該更新磁化率分布が算出される毎に平滑化するステップと、
を含む磁化率分布計算プログラム。
【国際調査報告】
2016076076.tif 000022
2016076076.tif 000023
2016076076.tif 000024
2016076076.tif 000025