JP2019162369A - 磁気共鳴撮影装置、感度分布算出方法、及び、感度分布算出プログラム - Google Patents

Abstract

【課題】被検体の形状や構造に依存しない高解像度コイル感度分布を高精度に算出する。【解決手段】ＭＲＩ装置は、複数のチャンネルからなる受信コイルを備え、前記受信コイルのチャンネル毎に被検体の核磁気共鳴信号を計測する計測部と、前記受信コイルのチャンネル毎の感度分布と、前記計測部が計測したチャンネル毎の核磁気共鳴信号から得たチャンネル画像とを用いて前記被検体の画像を作成する画像演算部と、を備え、前記画像演算部は、前記チャンネル画像と当該チャンネル画像を合成した合成画像とを用いて、チャンネル毎に、ｋ空間上で感度分布を算出する感度分布算出部を備える。【選択図】図９

Description

本発明は、核磁気共鳴を利用して検査対象の画像を取得する磁気共鳴撮影装置（以下、ＭＲＩ装置という）に関し、特に核磁気共鳴信号を受信する受信コイルの感度分布を算出する技術に関する。

ＭＲＩ装置では、静磁場内に置かれた被検体から核磁気共鳴によって発生する信号を受信コイルで受信し、その信号を処理することで被検体の画像を取得する。一般に受信コイルの感度は一様ではなく、感度分布が存在する。感度分布があると画像にシェーディングが生じるため、従来、感度分布の影響を補正により除去している。

一方、近年のＭＲＩでは、受信コイルによって感度分布が異なることを利用して、複数の感度分布の異なる受信コイルを用いて高速撮影する手法（パラレルイメージング）が広く採用されている。パラレルイメージングは、ｋ空間データを間引いて計測し、未計測点を、複数の受信コイルの感度分布を用いた演算によって復元し、折り返しのない画像を得る撮影方法であり、計測方法や信号復元方法が異なる種々の方法が提案されている。代表的なものとして、ｋ空間を間引いて計測し、計測空間（ｋ空間）上の演算によって、未計測のデータを推定して、データがフルに充填されたｋ空間データを作り、折り返しを発生させないようにする方法（ＳＭＡＳＨ法やＧＲＡＰＰＡ法など）、ｋ空間データから再構成した画像に生じた折り返しを、実空間上の演算で展開し、折り返しのない画像にする方法（ＳＥＮＳＥ法）がある（非特許文献１）。さらに複数のスライスを同時に励起して、得られた複数スライスからの信号を実空間上で信号分離する方法もある（ＳＭＳ法：Ｓｉｍｕｌｔａｎｅｏｕｓ Ｍｕｌｔｉｓｌｉｃｅ）。

上述した実空間上で信号分離する画像再構成では、演算に用いるコイル感度分布の精度が最終的に得られる画像の画質に大きく影響するため、コイル感度分布を高精度に求める必要がある。このため、被検体存在下でプリスキャンによって低解像度の感度分布算出用画像を取得し、この画像からコイル感度分布を算出し、さらに算出された低解像度の感度分布を補間拡大したり無信号／低信号領域の外挿を行ったりすることで、高解像度感度分布を得る手法が提案されている（特許文献１）。

Ｋｌａａｓ Ｐ．Ｐｒｕｅｓｓｍａｎｎ， ｅｔ ａｌ． Ｍａｇｎｅｔｉｃ Ｒｅｓｏｎａｎｃｅ ｉｎ Ｍｅｄｉｃｉｎｅ ４２：９５２−９６２（１９９９）

特開２００４−３４４３２７号公報

しかし上述した従来の感度分布の算出方法では、被検体の形状によって補間や外挿処理の精度が低下するという問題がある。それは、画像を構成する画素の輝度値は、被検体の形状や構造を反映したものであり、そのような画素の輝度値で構成される画像から算出した感度分布は、受信コイル自体の感度分布に、被検体の形状や構造に依存する要素が含まれることとなるためである。

そこで本発明は、被検体の形状や構造に依存しないコイル感度分布を高精度に算出すること、これによりコイル感度分布を用いたパラレルイメージング等の撮影画像の画質向上を図ることを課題とする。

上記課題を解決するため本発明では、ｋ空間上の演算により、コイル感度分布の周波数成分を求め、この周波数成分を用いて実空間データとしてのコイル感度分布を算出する。ｋ空間の周波数成分は、ｋ空間各位置の共通の周波数成分として求められるので、被検体の形状や構造に依存しない感度分布が得られる。

すなわち本発明のＭＲＩ装置は、複数のチャンネルからなる受信コイルを備え、前記受信コイルのチャンネル毎に被検体の核磁気共鳴信号を計測する計測部と、前記受信コイルのチャンネル毎の感度分布と、前記計測部が計測したチャンネル毎の核磁気共鳴信号から得たチャンネル画像とを用いて前記被検体の画像を作成する画像演算部と、を備え、前記画像演算部は、前記チャンネル画像と当該チャンネル画像を合成した合成画像とを用いて、チャンネル毎に、ｋ空間上で感度分布を算出する感度分布算出部を備える。

本発明によれば被検体の形状や構造に依存しない高精度の感度分布を得ることができる。それによりパラレルイメージング等のコイル感度分布を用いた撮影方法において、精度のよい画像を得ることができる。

ＭＲＩ装置の機能ブロック図。 実施形態のＭＲＩ装置の全体構成を示す図。 第一実施形態の撮影手順を示す図。 第一実施形態における感度分布算出用パルスシーケンスの一例を示す図。 感度分布算出の手順を示す図。 感度分布算出部の機能ブロック図。 図５のフローのマルチコイル合成ステップの詳細を示す図。 図５のフローの周波数成分算出ステップＳ１５３の詳細を示す図。 第一実施形態の処理を示す図。 第一実施形態により得られた高解像度感度分布を示す図。 第一及び第二実施形態とその変形例の処理の概要を示す図。 第三実施形態の撮影手順を示す図。

以下、図面を参照して、本発明のＭＲＩ装置の実施形態を説明する。

＜装置の実施形態＞
最初に後述する各実施形態に共通するＭＲＩ装置の実施形態を説明する。
本実施形態のＭＲＩ装置は、図１に示すように、核磁気共鳴信号を取得する計測部１０と、計測部が取得した核磁気共鳴信号（以下、ＭＲ信号という）を用いて画像再構成等の演算を行う演算部２０と、計測部１０及び演算部２０などＭＲＩ装置を構成する各部の動作を制御する制御部３０とを備える。制御部３０は、演算部２０の機能の一部であってもよい。さらにＭＲＩ装置は、ユーザーが計測部１０や制御部３０に撮影の開始／終了、撮影条件の設定の指令等を入力したり、演算部２０の処理結果をユーザーに提示したりするためのデバイスを備えたユーザーインターフェイス装置４０を備えることができる。

計測部１０は、一般的なＭＲＩ装置と同様の構成を有し、被検体が置かれる空間に均一な磁場を発生する静磁場発生部１１、静磁場に対し磁場勾配を与える傾斜磁場発生部１２、被検体の組織を構成する原子の原子核、典型的にはプロトン、を励起する高周波磁場パルス（以下、ＲＦパルスという）を照射する送信部１３、被検体から発生する核磁気共鳴信号を受信する受信部１４を備える。具体的には、図２に示すように、静磁場発生部１１として静磁場コイル１０２とその磁場均一性を補正するシムコイル１０４及びその電源部１１３、傾斜磁場発生部１２として傾斜磁場コイル１０３とその電源部１１３、送信部１３として送信コイル１０５と送信コイル１０５にＲＦ信号を送る送信機１０７、受信部として受信コイル１０６と受信コイル１０６が受信したＭＲ信号を検波する受信機１０８、を備えている。

静磁場発生部１１は、常伝導磁石方式または超伝導磁石方式の静磁場コイル１０２、或いは永久磁石を用いることができる。また磁場の向きによって垂直磁場方式、水平磁場方式、或いは磁場の向きが水平方向に対し傾斜した方式などがあり、そのいずれを採用してもよい。

傾斜磁場コイル１０３は、図２では一つのブロックで示しているが、互いに直交する３軸方向の傾斜磁場コイルからなり、それぞれがそれぞれの電源部１１２に接続されている。各軸の傾斜磁場の印加により、被検体の励起領域を選択したりＭＲ信号に位置情報を付与したりすることができる。

送信コイル１０５は、一般的には、シムコイル１０４及び傾斜磁場コイル１０３とともに静磁場空間内に固定されている。

受信コイル１０６は、図２では一つのブロックで示しているが、複数のチャンネルを持つ受信コイルで構成されており、チャンネル毎に受信を行う機能を持つ。複数チャンネルの受信コイルとしては、それぞれがチャンネルに対応するサブコイルを組み合わせた受信コイルや、複数のループコイルを所定の配置で配置した受信コイルなど公知の複数チャンネルを用いることができ、ここではサブコイル或いはループコイルがそれぞれ個々の受信コイルに相当する。

計測部１０は、制御部３０のもと、ユーザーが指定する撮像条件と所定のパルスシーケンスに従い、送信コイル１０５及び傾斜磁場コイル１０３を動作させて、静磁場内の被検体１０１にＲＦパルスを照射し、これにより被検体から発生するＭＲ信号を受信コイル１０６が受信する。制御部３０による制御のうち、パルスシーケンスによる計測部１０の制御は、シーケンス制御装置１１４により実現される。

演算部２０は、計測部１０が計測した核磁気共鳴信号を用いて、受信コイルの感度分布を算出する感度分布算出部２１と、感度分布算出部２１が算出した感度分布及び計測データであるｋ空間データを用いて画像を生成する画像生成部２３と、を備える。演算部２０の機能は、例えば、ＣＰＵ或いはＧＰＵとメモリを備えた計算機に、感度分布算出や画像生成の手順を定めたプログラムを搭載し実行させることによって実現することができる。但し、演算の一部或いは全部を、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）やＦＰＧＡ（Ｆｉｅｌｄ Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）等のハードウェアで構成することや、ＭＲＩ装置とは別の計算機やクラウド上の演算手段にて実現することも可能である。

制御部３０は、計測部１０を所定のパルスシーケンスに従って制御するシーケンス制御装置１１４を含み、計測部１０と演算部２０の動作を制御するものであり、その機能の一部或いは全部は、演算部２０が計算機のＣＰＵで実現される場合に、それと同じＣＰＵに制御プログラムを搭載して実現してもよいし、一部或いは全部を別のＣＰＵやハードウェアで実現することも可能である。

図２に示す構成では、ソフトウェアで実現する演算部２０及び制御部３０の機能は計算機１０９に組み込まれ、計算機１０９にはユーザーインターフェイス装置４０となるディスプレイ１１０や入力装置１１５、及び外部記憶装置１１１などが接続されている。

本実施形態のＭＲＩ装置は、演算部（感度分布算出部２１）２０が、ｋ空間上の演算により、コイル感度分布の周波数成分を求め、この周波数成分を用いて実空間データとしてのコイル感度分布を算出することが特徴であり、計測部１０が実行する撮影方法は特に断らない限り、公知のＭＲＩ技術を採用することができ、また画像生成部２３が画像を生成する手法や生成する画像の種類などは、撮影方法に応じた公知のものが採用できる。

以下、感度分布算出部２１の機能を中心に、感度分布算出方法の実施形態を説明する。必要に応じて、計測部１０及び画像生成部２３の機能にも言及する。

＜第一実施形態＞
本実施形態では、複数の受信コイルの感度分布の違いを利用して高速撮影を行うパラレルイメージングを例に、ＭＲＩ装置の動作を説明する。また一例として、三次元（３Ｄ）の感度分布用シーケンスで得たデータを用いて、二次元（２Ｄ）の感度分布を算出する場合を説明する。

パラレルイメージングの撮影は、図３に示すように、計測（Ｓ１００１）、感度分布算出（Ｓ１００２）、信号分離処理（Ｓ１００３）、及び画像表示（Ｓ１００４）の各ステップを含む。以下、各ステップを詳細に説明する。

［計測ステップＳ１００１］
計測ステップＳ１００１では、計測部１０が被検体画像を取得するためのパルスシーケンス（以下、本撮像シーケンスという）と感度分布用画像を取得するためのパルスシーケンス（以下、感度分布用シーケンスという）と、を含む撮像シーケンスを実行する。パルスシーケンスの種類は特に限定されないが、図４に感度分布用シーケンスの例として、典型的な３Ｄ−ＧｒＥ（グラディエントエコー）法のパルスシーケンス４００を示す。図中、上から順に、ＲＦは高周波磁場パルス４０２の照射タイミング及び核磁気共鳴信号４０７の受信タイミングを示し、Ｇｓ、Ｇｐ、Ｇｒは、それぞれ互いに直交する３軸方向の傾斜磁場パルスの極性／強度と印加タイミングを示す。三桁の数字は、同じ数字は同種のパルスであることを示し、三桁の数字の後にハイフンで続く数字は、同じ数字は同じ繰り返し時間内で印加されるパルスであることを示す。また、Ｇｓはスライス選択傾斜磁場、Ｇｐは位相エンコード傾斜磁場、Ｇｒは読み出し傾斜磁場であり、パルスの横に示す縦方向の矢印は、傾斜磁場強度がパルスシーケンスの繰り返しごとに変化することを示している。

図４の撮像シーケンスでは、まずＲＦパルス４０２とスライス選択傾斜磁場４０１を印加して、所定の領域を励起した後、スライス傾斜磁場Ｇｓのリフェイズパルス４０３と位相エンコードパルス４０４を印加する。このとき読み出し傾斜磁場Ｇｒのディフェイズパルス４０５を印加し、その極性を反転させたリフェイズパルスを印加しながら、核磁気共鳴信号４０７をサンプリングする。最後に位相エンコード方向にリフェイズパルス４０９を印加する。ＲＦパルス４０２印加からリフェイズパルス５０９印加までを所定の繰り返し時間ＴＲで繰り返し、画像再構成に必要な核磁気共鳴信号を収集する。サンプリングされた核磁気共鳴信号のデジタルデータがｋ空間データとなる。このシーケンスでは３Ｄ−データが得られる。

本撮像シーケンスと感度分布用シーケンスとは同じでも異なっていてもよいが、感度分布用シーケンスの場合には、スライスエンコード及び位相エンコードともにエンコードステップ数は少なく、低解像度の感度分布用画像を取得する。これにより、本撮像シーケンスとは別のプリスキャンとして感度分布用シーケンスを実行する場合に、プリスキャンに要する計測時間を短時間に抑えることができる。但し、低解像度の感度分布は、高解像度の本撮像の画像とは解像度や視野（ＦＯＶ）が一致しないため、コイル感度分布を画像再構成に適用するためには、後述する感度分布算出処理において、本撮像の画像に一致させる処理を行う。

［感度分布算出ステップＳ１００２］
感度分布算出ステップＳ１００２では、感度分布算出部２１が、計測部１０により収集された感度分布用のデータ（以下、Ｓｍａｐデータという）を用いて本撮像の画像（対象画像）に適用されるコイル感度分布を算出する。本実施形態はパラレルイメージングを行うので、複数の受信コイル（チャンネル）のそれぞれについてコイル感度分布を算出する。

感度分布算出ステップＳ１００２の詳細を、図５に示す。感度分布算出は、大きく分けて、各受信コイルの感度分布（画像空間データ）の合成Ｓ１５１、合成感度分布のｋ空間データへの変換Ｓ１５２、ｋ空間データを用いた各受信コイルの周波数成分算出Ｓ１５３、感度分布の周波数成分（ｋ空間）の実空間データへの変換Ｓ１５４、を含み、さらに必要に応じてマスク画像を乗算する処理Ｓ１５５を含む。これらの処理を実現するため感度分布算出部２１は、例えば、図６の機能ブロック図に示すように、スライス切り出し部２１１、マルチコイル画像合成部２１２、フーリエ変換部２１３、感度分布周波数成分算出部２１４、δ関数画像算出部２１５、畳込み積分部２１６、マスク画像乗算部２１７を備える。なお図６は感度分布算出部２１の一例に過ぎず、機能部の一部を適宜省略したり、他の機能部を追加したりすることも可能である。以下、各部の動作として、感度分布算出部２１が行う処理の詳細を説明する。

（合成処理Ｓ１５１）
合成処理Ｓ１５１の詳細を図７に示す。
合成処理Ｓ１５１では、まずフーリエ変換部２１３が３Ｄ−Ｓｍａｐデータをフーリエ変換し、ｋ空間データから画像データ（実空間データ）に変換した後（Ｓ１５１１）、スライス切り出し部２１１が、実空間データである３Ｄ−Ｓｍａｐ画像データから、２Ｄ画像データ（スライス）を切り出す（Ｓ１５１２）。このとき、切り出すスライスのＳｍａｐ画像データは、そのＦＯＶが対象画像のＦＯＶと同じＦＯＶとなるように切り出す。一般に感度分布用計測では、対象画像のＦＯＶより大きなＦＯＶで感度分布画像を取得するので、感度分布画像のＦＯＶを対象画像のＦＯＶに一致させておく。但し、画像のマトリクスサイズは、感度分布画像が対象画像よりも小さい低解像度画像である。

その後、マルチコイル画像合成部２１２が各チャンネルのＳｍａｐ画像データを合成し、１枚の合成Ｓｍａｐ画像を得る（Ｓ１５１４）。

なお各チャンネルのＳｍａｐ画像データを合成する際に、低解像度画像であるＳｍａｐ画像のＳＮＲ（信号ノイズ比）を向上するために、ノイズ除去処理を行ってもよい（Ｓ１５１３）。ノイズ除去処理としては、例えば、画素値の二乗和（Ｓｕｍ−ｏｆ−Ｓｑｕａｒｅ）を取って画素値とするＳｏＳ画像をガイド画像としてフィルタリングする手法（Ｇｕｉｄｅｄ Ｆｉｌｔｅｒｉｎｇ法）などを採用することができる。この手法は頭部などの低解像度画像において低ＳＮＲとなりやすい画像に効果的に適用することができ、精度の高い合成Ｓｍａｐ画像を得ることができる。但し、ノイズ除去処理は、上記フィルタリング手法に限らず公知の手法を採用でき、また省略することも可能である。

マルチコイル画像合成部２１２によるＳｍａｐ画像データの合成処理は、次式（１）で表すことができる。

式中、Ｓ_ＭＡＣ＝は合成Ｓｍａｐ画像、Ｃｉはｉ番目のチャンネルのコイル感度分布、Ｓｉはｉ番目のチャンネルのＳｍａｐ画像で、Ｃｉ及びＳｉの「ｉ」は、チャンネル数がｎの場合、１〜ｎのいずれかである。ここでＣｉは、マルチプルコイル合成のための仮に算出されるコイル感度マップであり、式（２）により算出することができる。

式中、Ｓ_ＳｏＳは、上述したフィルタリングのガイド画像に用いた二乗和ＳｏＳ画像と同じである。また「exp」のφ_ｒｅｆは参照チャンネルの位相画像であり、この項を入れることで参照チャンネルとの位相差分を位相情報として含む感度マップが得られる。参照チャンネルは、特に限定されないが、信号値が最大となるチャンネルであることが望ましい。

以上の各処理Ｓ１５１１〜Ｓ１５１４により、合成処理Ｓ１５１が完了する。

（フーリエ変換Ｓ１５２）
フーリエ変換部２１３は、マルチコイル画像合成部２１２が合成処理Ｓ１５１により得た合成Ｓｍａｐ画像をフーリエ変換し、周波数空間すなわちｋ空間のデータとする。

（周波数成分算出Ｓ１５３）
感度分布周波数成分算出部２１４は、合成Ｓｍａｐ画像（Ｓ_ＭＡＣ）のｋ空間データから感度分布の周波数成分を算出する。まず感度分布の周波数成分について説明する。

実空間において、実空間上の各位置のベクトルをｒとすると、各チャンネル画像Ｓｉ（ｒ）は、式（３）に示すように、各コイルの感度分布Ｃｉ（ｒ）と合成Ｓｍａｐ画像Ｓ_ＭＡＣ（ｒ）との積で表される。

式（３）の両辺をｋ空間データに変換すると、式（３）の関係は、式（４）の畳込み積分で表される。

式中、ｋはｋ空間位置ベクトルである（以下、同じ）。

式（４）で示すように、Ｓ_ｉ（ｋ）は、Ｓ_ＭＡＣ（ｋ）と感度分布の周波数成分ｃ_ｉ（ｋ）が畳込まれているものである。ここでＳ_ＭＡＣ（ｋ）は、実空間画像をｋ空間データに変換したものであって、各画素の輝度値の違い、すなわち被検体の形状や構造を周波数として表しているが、ｃ_ｉ（ｋ）は、ｋ空間の各位置に共通の周波数成分であり、被検体の形状や構造に依存しない。従って、式（４）を解くことで、被検体の形状や構造に依存しない感度分布の周波数成分ｃ_ｉ（ｋ）が求まる。

感度分布の周波数成分ｃ_ｉ（ｋ）を算出する処理は、具体的には、図８に示すように、受信コイルの周波数成分に相当する重み付け係数（カーネル）を算出するステップ（Ｓ１５３１）と、重み付け係数とδ関数とを畳み込み積分するステップ（Ｓ１５３４、Ｓ１５３５）と、畳み込み積分に用いるδ関数を算出するステップ（Ｓ１５３２）とを含む。また図８に示すフローは、感度分布算出の精度を高めるために、繰り返し演算（Ｓ１５３３）を行うとともに、被検体の存在しない領域を取り除くマスク処理ステップ（Ｓ１５３６）を追加している。なお、繰り返し演算は省略することも可能である。

まず重み付け係数の算出（Ｓ１５３１）について説明する。
上述した式（４）は、所定のカーネルサイズの重み係数を用いて、式（５）のように表すことができる。

式（５）において、ｋｘ、ｋｙはｋ空間の位置、ｂｘ、ｂｙはカーネル内のｋ空間位置、Ｂｘ、Ｂｙはカーネルサイズである。カーネルサイズは、特に限定されるものではないが、たとえば、５×５や９×９など小さいサイズとする。

式（５）は、行列表現すると式（６）となり、式（７）で示す擬似逆行列から、重みベクトルｗｉ（すなわち、重み係数）を算出することができる。

なおＳｍａｐ画像のサイズをそれぞれＮｘ、Ｎｙ（但しＮｘ及びＮｙは偶数）であるとすると、行列ＳＭＡＣのサイズは、カーネルサイズが奇数の場合、
（縦，横）＝｛（Ｎｘ−Ｂｘ＋１）×（Ｎｙ−Ｂｙ＋１）,Ｂｘ×Ｂｙ｝
となる。

畳込み積分部２１６は、こうして得られる重み係数ｗｉに、式（８）に示すように、所望のサイズのδ関数を畳込み積分し、式（４）に示す感度分布の周波成分ｃ_ｉ（ｋ）を得る（Ｓ１５３４またはＳ１５３５）。すなわち重み係数は、感度分布の周波成分であるがカーネルサイズが小さいものであるので、所望のサイズのδ関数を畳込み積分することで、ｋ空間での外挿を行って、所望サイズの周波数成分が得られる。

このため、δ関数画像算出部２１５は、所望のサイズですべての画素値を１とする画像を作成し、それをフーリエ変換することにより、所望サイズのδ関数を作成する（Ｓ１５３２）。例えば、低解像度の合成Ｓｍａｐ画像Ｓ_ＭＡＣのサイズが６４×６４の場合、サイズ６４×６４の画像（画素値＝１）をフーリエ変換することで６４×６４のδ関数を作成する。ここでは、スライス（２Ｄデータ）を対象としているので二次元のδ関数を作成する。図８に示すフローにおいて、繰り返し演算を省略する場合には、本撮像の高解像度画像のサイズに合わせたδ関数を作成し、図８で点線の矢印で示すように、これを式（８）に用いてもよい（Ｓ１５３５）。

繰り返し演算を行う場合（Ｓ１５３３−１の「ＹＥＳ」）には、ステップＳ１５３４で、低解像度のδ関数を用いて感度分布の周波数成分を算出した後、これを各チャンネルのＳｍａｐ画像を合成したときに式（１）において使用した感度マップ（式（２）によって求めた仮の感度マップ）に置き換えて、再度合成Ｓｍａｐ画像を作成し（Ｓ１５１４）、その後の処理（Ｓ１５２、Ｓ１５３）を繰り返す。処理を繰り返すことにより、最終的に得られる感度マップの精度を高めることができる。

繰り返しの終了（Ｓ１５３３−２）は、例えば、前回に得られた感度マップの二乗と今回に得られた感度マップの二乗との差や次式（９）で示す変化の割合などが所定の閾値以下になった時点とする、あるいは繰り返し回数を終了条件としてもよい。

式（９）において、ａは所定の閾値、例えば、「０．０１」、ｍはｍ回目の繰り返しであることを示す。

繰り返しが終了した後、本撮像の高解像度画像のサイズに合わせたδ関数を用いて、高解像度の感度分布の周波数成分を算出する（Ｓ１５３５）。

（実空間への変換：Ｓ１５４、マスク乗算：Ｓ１５５）
最後に、上記ステップＳ１５３５で得られた感度分布の周波数成分をフーリエ変換し、実空間データである感度分布（感度マップ）を得る。その後、感度分布の周辺のノイズを除去するために所定のマスク画像を乗算する。マスク画像乗算部２１７は、は、感度分布用画像として取得した低解像度画像からマスク画像を作成し、フーリエ変換後の感度分布画像に乗算する。なおマスク画像の作成は、３Ｄ感度画像からスライスを切り出した後（図７：１５１２）であれば、いつでも作成することができる。

以上説明したＳ１５１〜Ｓ１５４により、図３の感度分布算出処理Ｓ１００２が終了する。図９に本実施形態の処理を模式的に示す。またチャンネル数５の受信コイルを用いて実際にヒト頭部の画像から本実施形態によって算出した高解像度の感度分布画像を図１０に示す。図１０において、（ａ）は、各チャンネルの画像の画素値の二乗和（Ｓｕｍ−ｏｆ−Ｓｑｕａｒｅ）を取って画素値としたＳｏＳ画像、（ｂ）は一つのチャンネルの画像、（ｃ）は同チャンネルの感度分布（マスク乗算後）である。

［信号分離ステップＳ１００３、画像表示ステップＳ１００４］
画像生成部２３（図１）は、本撮像シーケンスの実行により、各受信コイルで取得した画像に対し、上述の感度分布算出処理で算出した各受信コイルの感度分布を用いて、公知のパラレルイメージングアルゴリズムを用いて、信号分離（折り返し除去）演算を行い、被検体画像を得る。本実施形態ではスライス毎に感度分布を得ているので、本撮像シーケンスで得た画像が３Ｄ画像であれば、対応するスライス位置の画像に対し、スライス毎に信号分離処理を行う。マルチスライス画像であれば、スライス切り出し（Ｓ１５１２）において、本撮像のスライスと同じスライス位置を切り出して感度分布算出を行い、対応するマルチスライス画像の各スライスに対し、対応するスライス位置の感度分布を用いて信号分離処理を行う。

演算部２０は、被検体画像を、例えば、ユーザーインターフェイス装置４０のディスプレイ１１０に表示させたり、通信等を通じて、ＭＲＩ装置から離れた場所にあるディスプレイに表示させたりすることができる。

以上、説明したように、本実施形態によれば、まず各受信コイルの感度分布画像及びそれらを合成した合成感度分布画像をｋ空間データとして、受信コイルの感度分布の周波成分を求め、この周波数成分に対し所望のサイズのδ関数を用いた畳込み演算を行うことにより、被検体の形状や構造の影響を排除した高解像度の感度分布画像を得ることができる。これにより受信コイルの感度分布を用いた実空間上の信号分離処理の精度を高め、精度のよいパラレルイメージングを行うことができる。

また本実施形態では、スライス毎の感度分布（２Ｄ−感度分布）を算出するので、感度分布算出の際に用いる最適なカーネルの設計が容易である。また本撮像の画像のうち所望のスライスのみを信号分離処理することも可能であり、また信号分離処理を簡易な演算で実現できる。

＜変形例１＞
第一実施形態では、本発明をｋ空間データを間引いて計測するパラレルイメージングに適用した例を説明したが、データを間引かずに計測し、複数のチャンネルの画像に感度分布を乗算して合成する撮影（ＭＡＣ合成）にも第一実施形態の感度分布算出手法は適用できる。その場合、図３における信号分離処理（Ｓ１００３）は、信号合成処理になる。
本変形例においても、高精度の高解像度感度分布を用いることで、合成画像の精度を向上することができる。

＜変形例２＞
第一実施形態では、計測部１０が感度分布用画像として、３Ｄ−Ｓｍａｐ画像を取得し、この画像データからスライス（２Ｄ画像データ）を切り出し、２Ｄの感度分布を算出したが、計測部１０がマルチスライス撮影により複数の２Ｄ−Ｓｍａｐ画像を取得し、それぞれの、或いは一部のスライスについて、２Ｄの感度分布を算出してもよい。第一実施形態と本変形例の処理の概略をまとめて図１１に示す。図では、第一実施形態の手法を点線で、変形例の手法を一点鎖線で囲って示している。

３Ｄ−Ｓｍａｐ画像からスライスを切り出す場合には、切り出すスライスの自由度が高いが、本変形例では、図６のスライス切り出し部２１１及び図５のステップＳ１５１２（各チャンネルの３Ｄ画像データからスライスを切り出すステップ）が不要となり、感度分布計算の負担を減らすことができる。

＜第二実施形態＞
第一実施形態では、計測した３Ｄ−感度分布画像データから、または２Ｄ−感度分布画像データから、スライス毎に２Ｄ−感度分布を算出したが、本実施形態は３次元の感度分布を一回の処理で取得することが特徴である。

３次元の感度分布を算出する処理は、感度分布用シーケンスで得た感度分布用画像が３Ｄ−データの場合とマルチスライスの２Ｄ−データの場合のいずれにも適用できるが、ここでは第一実施形態と同様に３Ｄ−データである場合を説明する。適宜、図５及び図７を参照する。

まず、３Ｄ−感度分布用シーケンス（例えば図４のパルスシーケンス）を実行し、３Ｄ−感度分布画像データ（Ｓｍａｐ画像）を得た後、この３Ｄ−感度分布画像データを合成し、３Ｄ−合成感度分布画像を得る（合成ステップＳ１５１）。本実施形態では、図７に示す合成ステップＳ１５１のうち、３Ｄデータの感度分布を取得するため、スライス切り出しステップ（図７：Ｓ１５１２）は不要である。また仮の感度分布を用いた合成は、スライス毎の感度分布ではなく、３Ｄ感度分布を用いる。

次いで３Ｄ−合成感度分布画像をフーリエ変換し（Ｓ１５２）、次式（１０）〜式（１２）より、ｋ空間における各受信コイルの重み係数ｗｉを算出する。式（１０）〜（１２）は、第一実施形態で説明した式（５）〜（７）に対応し、各受信コイルのｋ空間データが、重み係数と合成感度分布画像のｋ空間データとの畳込み積分になっていることを用いて、擬似逆行列演算によって、３Ｄの重み係数（カーネル関数）を求めるものである。

ここで、Ｓｍａｐ画像のサイズをそれぞれＮ_ｘ, Ｎ_ｙ, Ｎ_ｚ （偶数）とすると、行列Ｓ_ＭＡＣのサイズは、
（縦，横）＝ ( (N_x-B_x+1) × (N_y-B_y+1) × (N_z-B_z+1), B_x × B_y × B_z )
となる（カーネルサイズが奇数のとき）。

その後、この３Ｄ重み係数とδ関数との畳込み積分を行い、所望のサイズの感度分布の周波数成分を得る。δ関数は、それを作成するもとの画像が３Ｄ画像であって、もとの画像と同じサイズの３Ｄの関数である。この畳込み積分の結果をフーリエ変換することにより、実空間データである感度分布を得ることができる。この感度分布算出過程においても、最終的に得られる感度分布の精度向上のために、低解像度δ関数を用いて、カーネル算出と畳込み積分の処理を繰り返してもよく、繰り返しによって精度が高まった状態で、高解像度δ関数を用いて最後の畳込み積分を行い、感度分布を得てもよい。

本実施形態によれば、スライス毎に２Ｄの感度分布を求める場合（第一実施形態）に比べ、ステップＳ１５３１における最適なカーネルの設計の難易度は高くなるが、スライス間のずれが発生するおそれがなく、なめらかな感度分布を得ることができる。

なお本実施形態では、計測部１０が取得する感度分布用画像データが３Ｄ−データであるが、２Ｄ―データであっても、マルチスライスのデータを３Ｄ−データとして扱うことで同様に1回の処理で３Ｄ−感度分布を算出することができる（変形例１）。さらに３Ｄ−データから、本撮像に対応する所定の複数スライスを切り出して、それを３Ｄ−データとして扱い、１回の処理で３Ｄ−感度分布を算出してもよい（変形例２）。本実施形態とこれら変形例の処理を、第一実施形態における変形例の処理の概略を説明した図１１に併せて示す。図中、本実施形態とその変形例の処理はＬ型の矢印で示している。

＜第三実施形態＞
上述した第一及び第二実施形態では、感度分布用シーケンスを、本撮像シーケンスとは別に行ったが、感度分布算出は、本撮像シーケンスで得られる計測データの一部を用いて行うこと（セルフキャリブレーション）も可能である。

本実施形態の処理を図１２に示す。図１２において、図３と同じ処理は同じ符号で示し説明を省略する。本実施形態の処理は、本撮像の計測データを感度分布算出に利用すること以外は、第一または第二実施形態と同様であり、第一または第二実施形態で説明した変形例も同様に適用することができる。

本実施形態の処理では、計測Ｓ１０１０では、感度分布用シーケンスを省略し、本撮像シーケンスのみを実行する。その際、低周波領域のデータは間引かずに計測する。この間引かずに計測した低周波領域のデータを抽出し（Ｓ１０１１）、感度分布算出に用いる（Ｓ１００３）。低周波領域とは、例えばｋ空間の原点を含む低周波領域の数ライン程度の低解像度データであり、第一実施形態において、感度分布用シーケンスで取得したデータと同程度のデータである。その後、感度分布算出（Ｓ１００２）、信号分離（Ｓ１００３）、及び適宜画像表示することは第一実施形態と同様である。

本実施形態においても、図１１に示したように、感度分布算出には３Ｄ画像データ或いは２Ｄ画像データを用いて２Ｄの感度分布を算出してもよいし、３Ｄの感度分布を算出してもよい。本撮像のシーケンスに応じて適宜、いずれかの手法を採用すればよい。
本実施例によれば、セルフキャリブレーションを行うことにより、感度分布計測のための時間が不要になるため、体動の影響を軽減することができる。

以上、本発明の各実施形態を説明したが、本発明は周波数空間での感度分布演算を可能にする手段を備える装置であれば、これら実施形態に限定されるものではなく、図示した構造や工程に対し、要素を追加したり一部の要素を省略したりすることも可能である。

また本発明は、ｋ空間データを間引くパラレルイメージングだけでなく、複数スライス選択励起法（ＳＭＳ）パラレルイメージングや複数チャンネル画像の合成（ＭＡＣ合成）にも同様に適用することができる。一般的なＭＡＣ合成では、感度分布の二乗和を用いるため位相情報が消えてしまうが、本実施形態の手法で求めた感度分布を用いることで、位相情報を持つ画像（複素画像）を得ることができる。

１０・・・計測部、２０・・・演算部、２１・・・感度分布算出部、２３・・・画像生成部、３０・・・制御部、４０・・・ユーザーインターフェイス、２１１・・・スライス切り出し部、２１２・・・マルチコイル画像合成部、２１３・・・フーリエ変換部、２１４・・・感度分布周波数成分算出部、２１５・・・δ関数画像算出部、２１６・・・畳込み積分部、２１７・・・マスク画像乗算部。

Claims (12)

1. 複数のチャンネルからなる受信コイルを備え、前記受信コイルのチャンネル毎に被検体の核磁気共鳴信号を計測する計測部と、
   前記受信コイルのチャンネル毎の感度分布と、前記計測部が計測したチャンネル毎の核磁気共鳴信号から得たチャンネル画像とを用いて前記被検体の画像を作成する画像演算部と、を備え、
   前記画像演算部は、前記チャンネル画像と当該チャンネル画像を合成した合成画像とを用いて、チャンネル毎に、ｋ空間上で感度分布を算出する感度分布算出部を備えたことを特徴とする磁気共鳴撮影装置。
2. 請求項１に記載の磁気共鳴撮影装置であって、
   前記感度分布算出部は、前記チャンネル画像と前記合成画像とをｋ空間データに変換するｋ空間変換部と、
   前記ｋ空間データにおける、各チャンネルの重みを感度分布の周波数成分として算出する感度分布周波数成分算出部と、
   所定の撮像視野及び所定の分解能に対応するδ関数と、前記感度分布周波数成分とを畳込み積分する畳込み積分部と、
   前記畳込み積分後の周波数成分を実空間データに変換する実空間変換部と、を
   備えたことを特徴とする磁気共鳴撮影装置。
3. 請求項１に記載の磁気共鳴撮影装置であって、
   前記感度分布算出部は、前記計測部が計測したチャンネル毎の低解像度画像を用いて感度分布を算出することを特徴とする磁気共鳴撮影装置。
4. 請求項１に記載の磁気共鳴撮影装置であって、
   前記感度分布算出部は、前記計測部が計測したチャンネル毎の高解像度画像のうち低周波数領域のデータを用いて感度分布を算出することを特徴とする磁気共鳴撮影装置。
5. 請求項２に記載の磁気共鳴撮影装置であって、
   前記畳込み積分部は、前記δ関数として、前記計測部が計測した画像の撮像視野及び分解能に応じたδ関数を用いることを特徴とする磁気共鳴撮影装置。
6. 請求項１に記載の磁気共鳴撮影装置であって、
   前記感度分布算出部は、２次元スライス毎に感度分布を算出することを特徴とする磁気共鳴撮影装置。
7. 請求項１に記載の磁気共鳴撮影装置であって、
   前記感度分布算出部は、３次元の感度分布を一回の処理で算出することを特徴とする磁気共鳴撮影装置。
8. 受信コイルのチャンネル毎に計測した画像を用いてチャンネル毎の感度分布を算出する方法であって、
   前記チャンネル毎のチャンネル画像を合成して合成画像を作成するステップ、
   前記チャンネル画像と前記合成画像とをｋ空間データに変換するステップ、
   前記ｋ空間データにおける各チャンネルの感度分布周波数成分を算出するステップ、
   所定の撮像視野及び所定の解像度に対応するδ関数と前記感度分布周波数成分とを畳込み積分するステップ、及び
   前記畳込み積分後のデータを実空間データに変換するステップ、を含むことを特徴とする感度分布算出方法。
9. 請求項８に記載の感度分布算出方法であって、
   前記チャンネル画像は三次元画像であって、当該三次元画像からスライス画像を切り出すステップをさらに含み、スライス毎に感度分布を算出することを特徴とする感度分布算出方法。
10. 請求項８に記載の感度分布算出方法であって、
    前記畳込み積分するステップは、低解像度のδ関数を畳込み積分するステップと、高解像度のδ関数を畳込み積分するステップと、を含み、
    前記合成画像を作成するステップから前記畳込み積分するステップまでを、合成画像を作成するステップで用いる感度分布を更新しながら、複数回繰り返した後、高解像度のδ関数を畳込み積分するステップに進むことを特徴とする感度分布算出方法。
11. 請求項８に記載の感度分布算出方法であって、
    前記感度分布周波数成分を算出するステップは、前記感度分布周波数成分として、前記合成画像のｋ空間データから各チャンネルの重み係数を算出するステップを含むことを特徴とする感度分布算出方法。
12. コンピュータに、
    受信コイルのチャンネル毎のチャンネル画像を合成して合成画像を作成するステップ、
    前記チャンネル画像と前記合成画像とをｋ空間データに変換するステップ、
    前記ｋ空間データを用いて、各チャンネルの感度分布周波数成分を算出するステップ、
    所定の撮像視野及び所定の解像度に対応するδ関数を設定するステップ、
    前記感度分布周波数成分と前記δ関数とを畳込み積分するステップ、及び
    前記畳込み積分後のデータを実空間データに変換するステップ、
    を実行させる感度分布算出プログラム。

Images