(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2024-04-04
(45)【発行日】2024-04-12
(54)【発明の名称】磁気共鳴画像高速再構成法及び磁気共鳴イメージング装置
(51)【国際特許分類】
A61B 5/055 20060101AFI20240405BHJP
【FI】
A61B5/055 376
(21)【出願番号】P 2020025707
(22)【出願日】2020-02-18
【審査請求日】2022-12-23
(31)【優先権主張番号】P 2019040207
(32)【優先日】2019-03-06
(33)【優先権主張国・地域又は機関】JP
【新規性喪失の例外の表示】特許法第30条第2項適用 ウェブサイト:https://www.ieice.org/ken/paper/2019030661MB/ 公開日:平成31年2月26日 集会名:電子情報通信学会研究会(第36回IBISML研究会) 開催日:平成31年3月6日 授業名:国立大学法人九州大学システム情報科学研究院情報学専攻修士2年次の演習授業「情報学講究」 公開日:令和1年7月2日 ウェブサイト:https://link.springer.com/article/10.1007/s10334-019-00755-1 https://link.springer.com/journal/10334/32/1/suppl(Pages235-371の分) 公開日:2019年9月24日 ウェブサイト:https://esmrmb2019.smart-abstract.com/eposter/#/posters/124 公開日:2019年10月3日 集会名:European Society for Magnetic Resonance in Medicine and Biology 2019 Congress(ESMRMB 2019 Congress) 開催日:2019年10月5日
(73)【特許権者】
【識別番号】504145342
【氏名又は名称】国立大学法人九州大学
(73)【特許権者】
【識別番号】304024865
【氏名又は名称】学校法人杏林学園
(74)【代理人】
【識別番号】100095407
【氏名又は名称】木村 満
(74)【代理人】
【識別番号】100168114
【氏名又は名称】山中 生太
(74)【代理人】
【識別番号】100133592
【氏名又は名称】山口 浩一
(74)【代理人】
【識別番号】100162259
【氏名又は名称】末富 孝典
(72)【発明者】
【氏名】竹内 純一
(72)【発明者】
【氏名】實松 豊
(72)【発明者】
【氏名】川喜田 雅則
(72)【発明者】
【氏名】北崎 自然
(72)【発明者】
【氏名】久原 重英
【審査官】永田 浩司
(56)【参考文献】
【文献】米国特許出願公開第2017/0213321(US,A1)
【文献】Chang Min Hyun, et al.,Deep learning for undersampled MRI reconstruction,Phys. Med. Biol.,2018年,63,1-15
【文献】Seungjun Nah, et al.,Deep Multi-scale Convolutional Neural Network for Dynamic Scene Deblurring,CVPR,2017年,pp.3883-3891
【文献】Jasper Schoormans, et al.,Reconstruction of sparsely sampled Magnetic Resonance Imaging measurements with a convolutional neural network,1st Conference on Medical Imaging with Deep Learning,2018年,1-3
(58)【調査した分野】(Int.Cl.,DB名)
A61B 5/055
(57)【特許請求の範囲】
【請求項1】
磁気共鳴により観測される原信号
から、低周波数領域のサンプリング点が高周波数領域のサンプリング点よりも密となるように、サンプリング点を間引きして第1のk空間データを取得するデータ取得ステップと、
多重解像度畳み込みニューラルネットワークを用いて、前記第1のk空間データから高解像度画像を復元するデータ復元ステップと、を含み、
前記多重解像度畳み込みニューラルネットワークでは、
前記第1のk空間データから
、前記第1のk空間データよりも
サイズが小さく
かつ高周波成分が除去された第2のk空間データを生成し、
前記第1のk空間データ及び前記第2のk空間データ各々から得られる画像を、超解像畳み込みニューラルネットワークを用いて超解像化し、
超解像化された画像を、サイズを合わせて重ね合わせることにより、前記高解像度画像を復元する、
磁気共鳴画像高速再構成法。
【請求項2】
前記データ復元ステップに先立って、
高精細なk空間データを教師データとし、前記教師データが間引きされて生成されたk空間データを前記第1のk空間データとして、前記多重解像度畳み込みニューラルネットワークを学習させる深層学習ステップを含む、
請求項1に記載の磁気共鳴画像高速再構成法。
【請求項3】
前記深層学習ステップでは、
所定の空間周波数よりも低周波成分のデータ全体を残し、前記空間周波数よりも高周波成分についてはランダムにデータを残すアンダーサンプリングフィルタを用いて、前記教師データを間引きすることにより、学習用の前記第1のk空間データを生成する、
請求項2に記載の磁気共鳴画像高速再構成法。
【請求項4】
前記深層学習ステップでは、
前記教師データと、前記多重解像度畳み込みニューラルネットワークを用いて復元された前記高解像度画像に対応する第3のk空間データとの誤差に基づく誤差逆伝播法を用いて学習を行う、
請求項2又は3に記載の磁気共鳴画像高速再構成法。
【請求項5】
磁気共鳴により観測される原信号
から、低周波数領域のサンプリング点が高周波数領域のサンプリング点よりも密となるように、サンプリング点を間引きして第1のk空間データを取得するデータ取得部と、
多重解像度畳み込みニューラルネットワークを用いて、前記第1のk空間データから高解像度画像を復元するデータ復元部と、
を備え、
前記多重解像度畳み込みニューラルネットワークは、
前記第1のk空間データから
、前記第1のk空間データよりも
サイズが小さく
かつ高周波成分が除去された第2のk空間データを生成し、
前記第1のk空間データ及び前記第2のk空間データ各々から得られる画像を、超解像畳み込みニューラルネットワークを用いて超解像化し、
超解像化された画像を、サイズを合わせて重ね合わせることにより、前記高解像度画像を復元するように構成されている、
磁気共鳴イメージング装置。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、磁気共鳴画像高速再構成法及び磁気共鳴イメージング装置に関する。
【背景技術】
【0002】
磁気共鳴画像法(Magnetic Resonance Imaging; MRI)は生体内の情報を得る有力な手段の1つである。しかし、MRI装置の検査には、数十分から1時間ほどの時間が必要となる。そこで、患者の負担軽減のため、撮像時間のさらなる短縮が求められている。この10年間、情報処理によるMRIの高速化と高精細化に貢献してきたのが圧縮センシング(Compressed Sensing; CS)である(例えば、特許文献1、非特許文献1参照)。圧縮センシングは、アンダーサンプリングを行って撮像時間を削減しつつ、計算により少ない観測データから元の画像を復元し再構成する技術である。
【先行技術文献】
【特許文献】
【0003】
【非特許文献】
【0004】
【文献】J. Yang, J. Wright, T. Huang, and Y. Ma, "Image super-resolution via sparse representation," IEEE Trans. Image Processing, vol. 19, no. 11, pp. 2861-2873, November 2010.
【発明の概要】
【発明が解決しようとする課題】
【0005】
しかし、圧縮センシングのアルゴリズムは観測数の3乗オーダーの計算量を必要とする。このため、観測データを少なくして削減した撮像時間よりも、画像を復元し再構成する再構成時間の方が超えてしまうといった不都合が生ずる場合がある。
【0006】
本発明は、上記実情の下になされたものであり、撮像時間及び再構成時間を双方短縮しつつ、高解像度画像を復元することができる磁気共鳴画像高速再構成法及び磁気共鳴イメージング装置を提供することを目的とする。
【課題を解決するための手段】
【0007】
上記目的を達成するために、本発明の第1の観点に係る磁気共鳴画像高速再構成法は、
磁気共鳴により観測される原信号から、低周波数領域のサンプリング点が高周波数領域のサンプリング点よりも密となるように、サンプリング点を間引きして第1のk空間データを取得するデータ取得ステップと、
多重解像度畳み込みニューラルネットワークを用いて、前記第1のk空間データから高解像度画像を復元するデータ復元ステップと、を含み、
前記多重解像度畳み込みニューラルネットワークでは、
前記第1のk空間データから、前記第1のk空間データよりもサイズが小さくかつ高周波成分が除去された第2のk空間データを生成し、
前記第1のk空間データ及び前記第2のk空間データ各々から得られる画像を、超解像畳み込みニューラルネットワークを用いて超解像化し、
超解像化された画像を、サイズを合わせて重ね合わせることにより、前記高解像度画像を復元する。
【0008】
この場合、前記データ復元ステップに先立って、
高精細なk空間データを教師データとし、前記教師データが間引きされて生成されたk空間データを前記第1のk空間データとして、前記多重解像度畳み込みニューラルネットワークを学習させる深層学習ステップを含む、
こととしてもよい。
【0009】
また、前記深層学習ステップでは、
所定の空間周波数よりも低周波成分のデータ全体を残し、前記空間周波数よりも高周波成分についてはランダムにデータを残すアンダーサンプリングフィルタを用いて、前記教師データを間引きすることにより、学習用の前記第1のk空間データを生成する、
こととしてもよい。
【0010】
前記深層学習ステップでは、
前記教師データと、前記多重解像度畳み込みニューラルネットワークを用いて復元された前記高解像度画像に対応する第3のk空間データとの誤差に基づく誤差逆伝播法を用いて学習を行う、
こととしてもよい。
【0011】
本発明の第2の観点に係る磁気共鳴イメージング装置は、
磁気共鳴により観測される原信号から、低周波数領域のサンプリング点が高周波数領域のサンプリング点よりも密となるように、サンプリング点を間引きして第1のk空間データを取得するデータ取得部と、
多重解像度畳み込みニューラルネットワークを用いて、前記第1のk空間データから高解像度画像を復元するデータ復元部と、
を備え、
前記多重解像度畳み込みニューラルネットワークは、
前記第1のk空間データから、前記第1のk空間データよりもサイズが小さくかつ高周波成分が除去された第2のk空間データを生成し、
前記第1のk空間データ及び前記第2のk空間データ各々から得られる画像を、超解像畳み込みニューラルネットワークを用いて超解像化し、
超解像化された画像を、サイズを合わせて重ね合わせることにより、前記高解像度画像を復元するように構成されている。
【発明の効果】
【0012】
本発明によれば、処理時間が圧縮センシングよりも短い多重解像度畳み込みニューラルネットワークを用いて観測点数の少ない第1のk空間データをスケールダウンし、多重化して超解像処理を行って得られる画像を用いて画像の再構成を行う。このため、観測点数を少なくして撮像時間を短縮しつつ、多重解像度畳込みニューラルネットワークを用いて再構成時間を短縮して、高周波成分だけでなく低周波成分も損なうことなく画像を復元することができる。この結果、撮像時間及び再構成時間を双方短縮しつつ、高解像度画像を復元することができる。
【図面の簡単な説明】
【0013】
【
図1】本発明の実施の形態に係る磁気共鳴イメージング装置の構成を示すブロック図である。
【
図2】(A)は、被検体に印加されるXY方向の傾斜磁場を示す図である。(B)は、Z軸方向の傾斜磁場を示す図である。
【
図3】(A)は、2次元のk空間データを示す図である。(B)は、スライス画像を示す図である。
【
図4】多重解像度畳み込みニューラルネットワークの構成を示すブロック図である。
【
図5】超解像畳み込みニューラルネットワークの構成を示す模式図である。
【
図7】アンダーサンプリングフィルタの一例を示す図である。
【
図8】誤差逆伝搬法を用いた深層学習によるパラメータの更新の様子を示す模式図である。
【
図9】データ復元部のハードウエアの構成を示すブロック図である。
【
図10】本発明の実施の形態に係る磁気共鳴画像高速再構成法を示すフローチャートである。
【
図11】データ取得部で取得されたk空間データから再構成された断層画像(T1強調画像)を示す図である。
【
図12】データ復元部で復元されたk空間データから再構成された断層画像(T1強調画像)を示す図である。
【
図13】PSNRを測定する低周波領域を示す模式図である。
【
図14】PSNRを測定する中間領域を示す模式図である。
【
図15】PSNRを測定する高周波領域を示す模式図である。
【
図16】T1強調画像の復元の評価結果を示す図である。
【
図17】(A)は、頭頂方向から見たときに得られるMIP画像の元画像の一例である。(B)は、頭頂方向から見たときに得られるMIP画像の多重解像度畳み込みニューラルネットワークによる復元画像である。(C)は、頭頂方向から見たときに得られるMIP画像をアンダーサンプリングフィルタにより劣化させた劣化画像である。
【
図18】(A)は、
図17(A)の画像の白枠部分である。(B)は、
図17(B)の画像の白枠部分である。(C)は、
図17(C)の画像の白枠部分である。
【
図19】(A)は、横方向から見たときに得られるMIP画像の元画像の一例である。(B)は、横方向から見たときに得られるMIP画像の多重解像度畳み込みニューラルネットワークによる復元画像である。(C)は、横方向から見たときに得られるMIP画像をアンダーサンプリングフィルタにより劣化させた劣化画像である。
【
図20】(A)は、
図19(A)の画像の白枠部分である。(B)は、
図19(B)の画像の白枠部分である。(C)は、
図19(C)の画像の白枠部分である。
【
図21】(A)は、正面方向から見たときに得られるMIP画像の元画像の一例である。(B)は、正面方向から見たときに得られるMIP画像の多重解像度畳み込みニューラルネットワークによる復元画像である。(C)は、正面方向から見たときに得られるMIP画像をアンダーサンプリングフィルタにより劣化させた劣化画像である。
【
図22】(A)は、
図21(A)の画像の白枠部分である。(B)は、
図21(B)の画像の白枠部分である。(C)は、
図21(C)の画像の白枠部分である。
【
図23】劣化画像から復元画像への復元度を比較する表である。
【発明を実施するための形態】
【0014】
以下、本発明の実施の形態について図面を参照して詳細に説明する。各図面においては、同一又は同等の部分に同一の符号を付す。
【0015】
図1に示すように、磁気共鳴イメージング装置1は、データ取得部10と、データ復元部20と、を備える。
【0016】
データ取得部10は、静磁場の中に置かれた被検体(生体)に一連のシーケンスに従って高周波磁場及び傾斜磁場を印加して、被検体から発生する核磁気共鳴信号(原信号)を検出する。核磁気共鳴信号は、被検体内の水素密度分布を示す空間周波数の領域(k空間)の信号として観測される。データ取得部10は、磁気共鳴により観測される原信号に基づく第1のk空間データとしてのk空間データ15を取得する。
【0017】
具体的には、データ取得部10は、
図2(A)に示すように、被検体を+Z方向の静磁場B
0に置いた状態で、静磁場に垂直な方向から90°パルス及び180°パルスなどのRF(Radio Frequency)パルスを被検体に印加する。これにより、被検体の水素原子核が励起されて、核磁気共鳴信号(エコー信号)を取り出すことができる。データ取得部10は、この核磁気共鳴信号(エコー信号)をサンプリングして取得する。データ取得部10は、k空間データ15を、サンプリング点を間引きつつ取得することができる。
【0018】
データ取得部10において、RFパルスは、
図2(B)に示すように、Z軸方向に勾配を有するスライス傾斜磁場B(Z)と同時に被検体に印加される。通常は、Z軸方向が被検体の体軸方向となる。スライス位置Zに対応する静磁場強度をB(Z)とすると、スライス位置Z
0を選択励起するため、周波数γB(Z
0)を持つ正弦波のRFパルスを印加する。ここで、γは磁気回転比である。
【0019】
さらに、
図2(A)に示すように、データ取得部10において、被検体には、Y軸方向に傾きGyを有する位相エンコード傾斜磁場と、X軸方向に傾きGxを有する周波数エンコード傾斜磁場とが加えられる。さらに、核磁気共鳴信号(エコー信号)を観測するタイミングで、X軸方向に傾きGxを有する周波数エンコード傾斜磁場が再度印加される。上述のようなシーケンスが位相エンコード傾斜磁場の強度(傾きGy)を変えながら繰り返され、その度に核磁気共鳴信号(エコー信号)が観測される。
【0020】
上記シーケンスにおいて、観測される複数の核磁気共鳴信号(エコー信号)により、2次元のk空間データが生成される。
図3(A)に示すように、観測される核磁気共鳴信号(エコー信号)は、k空間では、点線で示される空間周波数Kyが同じ直線の領域に対応するデータとなる。
【0021】
k空間における位置座標(Kx,Ky)は、周波数エンコード傾斜磁場の傾きGx、位相エンコード傾斜磁場の傾きGyで表現すれば(γGxt,γGyt)となる。ここで、tは、位相エンコード傾斜磁場の傾きGy、周波数エンコード傾斜磁場の傾きGxの印加時間である。核磁気共鳴信号(エコー信号)は、位相エンコード傾斜磁場の傾きGyの強度を変えながら観測されるため、それぞれの信号はKy軸の位置が異なるデータとなる。このようにして、スライス位置Z0における2次元のk空間データが構成される。
【0022】
データ取得部10は、このような計測を、RFパルスの周波数を変更しながら、すなわちスライス位置Zを変更しながら繰り返し行う。これにより、スライス位置Zが異なる2次元k空間データが観測される。なお、このk空間データを逆フーリエ変換すれば、そのZ位置におけるXY面の2次元画像が得られる。その画像データを
図3(B)に示すように並べれば、被検体の3次元の画像データが得られるようになる。
【0023】
なお、データ取得部10でデータが取得される時間、すなわち撮像時間を短縮する場合には、シーケンスの繰り返し数又は観測するスライス位置Zの数を少なくする必要がある。シーケンスの繰り返し数を少なくすれば、
図3(A)に点線で示されるエコー信号の本数が少なくなり、スライス位置Zの数を少なくすれば、
図3(B)に示される2次元画像の枚数が少なくなる。
【0024】
図1に戻り、データ復元部20は、データ取得部10で取得されたk空間データ15を入力する。データ復元部20は、多重解像度畳み込みニューラルネットワークを用いて、k空間データ15から高解像度画像25を復元する。また、畳み込みニューラルネットワークは、もともと人間の脳の視覚に関わる構造を参考にして考案され、名前の通り畳み込み層を持つことを特徴とする。
【0025】
データ復元部20は、多重解像度畳み込みニューラルネットワーク21を用いてk空間データ15を高解像化し、高解像度画像25を生成して出力する。多重解像度畳み込みニューラルネットワーク21は、多重解像度(Multi-Resolution)を、畳み込みニューラルネットワーク(Convolutional Neural Network; CNN)で実現していることから、MRCNN21ともいう。
【0026】
MRCNN21は、入力したk空間データ15から、サイズの異なる(スケールダウンした)k空間データ(第2のk空間データ)を複数生成する。MRCNN21は、
図4に示すように、まず、k空間データ15を1/2、1/4にリサイズして低周波成分を取得し、2つのk空間データ(1/2データ、1/4データ)を生成する。k空間データ15、1/2データ及び1/4データはそれぞれ逆フーリエ変換される。例えば、入力されたk空間データ15が逆フーリエ変換されて生成される入力画像が、256×256
ピクセルのデータであった場合には、128×128
ピクセルの画像と、64×64
ピクセルの画像とが生成される
。
【0027】
続いて、サイズの異なる画像データそれぞれに対して超解像畳み込みニューラルネットワーク(SRCNN)35を用いた画像復元を行う。超解像とは、与えられた低解像度画像からより解像度の高い高解像度画像を得る手法のことである。ここでは、
図4に示すように、入力画像、1/2データに対応する画像、1/4データに対応する画像に対してそれぞれSRCNN35により超解像処理が行われる。
【0028】
SRCNN35は、
図5に示す3層の畳み込み層で構成され、大きく分けて3つのステップを含んでいる。第1のステップは、パッチ分解・抽出であり、第2のステップは、辞書学習・パッチ推定(非線形写像)であり、第3のステップは、データの再構成である。
【0029】
(第1のステップ)
第1のステップでは、パッチ分解・特徴量の抽出などが行われる。第1のステップでの写像(特徴マップ)を以下の関数F
1で表す。パッチとは画像を細かく分解したものである。パッチ分解することで2次元画像データのサイズが小さくなる。
【数1】
ここで、W
1は畳み込みのためのフィルタ群(重みフィルタ)であり、B
1は、バイアスである。W
1のサイズをc×f
1×f
1×n
1とする。cは入力データのチャネル数であり、f
1は、フィルタの空間的なサイズ(縦×横)を表し、例えば9である。n
1は、特徴マップの数であり、例えば64である。
【0030】
(第2のステップ)
第2のステップでは、低解像度パッチを高い解像度パッチに非線形的に写像し、圧縮センシングの計算を行う。データのサイズが小さいため、このときの圧縮センシングについては、短時間のうちに行うことができる。第2のステップである非線形写像に相当する写像(特徴マップ)を関数F
2で表す。
【数2】
ここで、W
2は重みフィルタを表す。W
2のサイズは、n
1×1×1×n
2である。n
2は、例えば32である。B
2は、バイアスを表すn
2次元ベクトルである。
【0031】
第3のステップでは、高解像度パッチベクトルを結合し、推定高解像度データdとして出力する。第3のステップである写像は、以下の関数Fで表すことができる。
【数3】
ここで、W
3は重みフィルタを表す。W
3のサイズは、n
2×f
3×f
3×cである。f
3は、例えば5である。B
3は、バイアスを表すc次元ベクトルである。ここで出力されるデータが推定高解像度データdとなる。
【0032】
ここで、
図4に示すように、k空間データ15に対応する入力画像が入力されるSRCNN35で推定される推定高解像度データdをd1とする。また、1/2データに対応する画像が入力されるSRCNN35で推定される推定高解像度データdをd2とする。さらに、1/4データに対応する画像が入力されるSRCNN35で推定される推定高解像度データdをd3とする。SRCNN35で復元された推定高解像度データd3は、2つの畳み込み層(CONV)を経て、アップサンプリングされ、推定高解像度データd2と重ね合わせられる。ここで、アップサンプリングとは、推定高解像度データを拡大する処理であり、アンプーリング(プーリングの逆の処理)とも呼ばれる。例えば、重ね合わせられた推定高解像度データd2+d3は、2つの畳み込み層(CONV)を経てアップサンプリングされ、k空間データ15に対応する画像から生成された推定高解像度データd1と重ね合わせられる。この重ね合わされた推定高解像度データd1+d2+d3が、2つの畳み込み層(CONV)を経て、復元された高解像度画像として復元される。この高解像度画像25は、学習時にフーリエ変換され、k空間データとして出力される。なお、各段の2つの畳み込み層のフィルタサイズは、3×3×2×2又は3×3×2×1である。
【0033】
データ復元部20において、MRCNN21によりk空間データ15から得られる画像の高解像化を行うためには、MRCNN21のパラメータを最適化する必要がある。このため、データ復元部20は、深層学習により、MRCNN21を学習させ、そのパラメータを最適化する学習部22を備えている。
【0034】
学習部22は、教師データ50を有している。教師データ50としては、間引きすることなくサンプリングされた被検体の高精細な断層画像をフーリエ変換して得られるk空間データが用いられる。学習部22は、教師データ50の一部が間引きされて生成されたデータを、k空間データ15としてMRCNN21に入力する。
【0035】
学習部22は、アンダーサンプリングフィルタ40を有している。アンダーサンプリングフィルタ40は、データからサンプルを間引いて解像度を下げるフィルタである。
図7では、k空間において、白い部分がデータを取得する部分であり、黒い部分がデータを間引く部分である。
図7に示すように、アンダーサンプリングフィルタ40は、k空間データにおける中間の空間周波数よりも低周波成分(白い領域)のデータ全体を残しつつ、空間周波数よりも高周波成分についてはランダムにデータを残すような形状のフィルタとなっている。アンダーサンプリングフィルタ40は、全体のデータ取得率を10%としている。例えば、低周波成分と高周波成分とのデータ取得率を1:1とすることができるが、この割合に特に限定されるものではない。
【0036】
このアンダーサンプリングフィルタ40は、低周波成分にエネルギーのほとんどが中心に集中しており、画像全体を高精度に復元することができるように構成されている。また、アンダーサンプリングフィルタ40は、また、高周波成分は、病変を発見する上で重要であるため、ポアソンディスクサンプリングによりランダムに高周波成分を取得するように構成されている。
【0037】
学習部22は、アンダーサンプリングフィルタ40を通過して間引きされたk空間データ15をMRCNN21に入力する。MRCNN21は、高解像度画像25を出力する。学習部22は、この場合に、MRCNN21で生成される高解像度画像25がフーリエ変換部41でフーリエ変換されて生成されるk空間データが教師データ50に近づくようにMRCNN21のパラメータ調整を行う。
【0038】
学習部22は、パラメータ設定部45を備えている。パラメータ設定部45は、教師データ50と、教師データ50が間引きされて生成されたk空間データ15から復元された高解像度画像25に対応するk空間データ(第3のk空間データ)との誤差に基づいて、誤差逆伝搬法を用いて深層学習を行う。誤差逆伝搬法は、入力層、中間層、出力層からなるネットワークに対し、出力層から入力層にかけて誤差の勾配を逆伝搬させることで各層の重みフィルタW1,W2,W3とバイアスB1,B2,B3といったパラメータを更新する教師付き学習アルゴリズムである。
【0039】
各層の重みについてパラメータの更新の様子を説明する。出力層付近に着目したネットワークの構造を
図8に示す。N層第kニューロンへの入力、N層第kニューロンの出力をそれぞれx
k
(N),y
k
(N)とする。また、N-1層第jニューロンからN層第kニューロンへの重みをw
jk
(N)とする。N層第kニューロンの活性化関数をg
k
(N)とおくと、x
k
(N),y
k
(N)は、ニューロンモデルから次の式で表すことができる。
【数4】
【数5】
【0040】
あるネットワークの入力に対する出力層の出力(高解像度画像25から得られるk空間データ)と教師データ50との誤差を表す関数を誤差関数L(y
1
(N),・・・,y
k
(N),・・・,y
K
(N))とする。まず、誤差関数LのN-1層第jニューロンからN層第kニューロンへの重みw
jk
(N)に関する勾配は、式(4)、(5)と連鎖律によって式(6)で表される。
【数6】
式(4)より、以下の式が成立する。
【数7】
【0041】
ただし、
【数8】
のように定義している。このように、誤差関数Lの勾配は、式(7)に示すように、出力側の層による項(式(8))と、出力層の一層前の出力y
j
(N-1)から計算することができる。
【0042】
次に、N-1層(中間層)の重みに関する勾配を考える。同様に、誤差関数LのN-2層第iニューロンからN-1層第jニューロンへの重みw
ij
(N-1)に関する勾配は、次式で表される。
【数9】
ここで、
【数10】
は、出力層から逆伝播された
【数11】
から計算することができる。このように誤差関数を入力層に向かって逆伝播させていき、各層のパラメータに関する誤差関数の勾配を求めていく方法が、誤差逆伝搬法である。
【0043】
パラメータの更新は、最急降下法の考え方に基づいて行う。最急降下法では誤差関数Lを用いて次式のように各パラメータを更新する。ただし、tは、反復回数を表す。
【数12】
重みの更新量は、誤差関数Lの勾配に学習係数ηをかけたものを重みの更新量とする。
【0044】
最急降下法では学習率が低いと収束が遅くなってしまう。また、tに対して勾配が振動するような場合も収束が遅くなる。これを解決するためにモーメンタム法と呼ばれる手法を用いるようにしてもよい。モーメンタム法における更新式を以下に示す。
【数13】
【数14】
μをモーメンタムと呼ぶ。上式(12)は、一回前の反復時の情報を用いて勾配の振動を抑制している。また勾配が振動せず同じ方向の場合は更新量を増加させる役割を持つ。なお、勾配法としてはadam法を用いるようにしてもよい。
【0045】
図9に示すように、データ復元部20は、制御部61、主記憶部62、外部記憶部63、操作部64をハードウエア構成として備えている。主記憶部62、外部記憶部63、操作部64はいずれも内部バス60を介して制御部61に接続されている。
【0046】
制御部61は、CPU(Central Processing Unit)等から構成されている。このCPUが、外部記憶部63に記憶されているプログラム69に従ってMRCNN21及び学習部22の処理を実行することにより、
図1に示すデータ復元部20が実現される。
【0047】
主記憶部62は、RAM(Random-Access Memory)等から構成されている。主記憶部62には、外部記憶部63に記憶されているプログラム69がロードされる。この他、主記憶部62は、制御部61の作業領域(データの一時記憶領域)として用いられる。
【0048】
外部記憶部63は、フラッシュメモリ、ハードディスク、DVD-RAM(Digital Versatile Disc Random-Access Memory)、DVD-RW(Digital Versatile Disc ReWritable)等の不揮発性メモリから構成される。外部記憶部63には、制御部61に実行させるためのプログラム69があらかじめ記憶されている。また、外部記憶部63は、制御部61の指示に従って、このプログラム69の実行の際に用いられるデータを制御部61に供給し、制御部61から供給されたデータを記憶する。
【0049】
操作部64は、キーボード及びマウスなどのポインティングデバイス等と、キーボードおよびポインティングデバイス等を内部バス60に接続するインターフェイス装置から構成されている。操作部64を介して、操作者が操作した内容に関する情報が制御部61に入力される。
【0050】
表示部65は、CRT(Cathode Ray Tube)またはLCD(Liquid Crystal Display)などから構成され、操作者が操作情報を入力する場合は、操作用の画面が表示される。
【0051】
入出力部66は、データ取得部10からk空間データ15を入力し、高解像度画像25を出力するインターフェイスとなる。
【0052】
データ復元部20では、制御部61がプログラム69を実行することにより、各種機能を実現する。
【0053】
なお、データ取得部10において演算を行う部分のハードウエア構成も、
図9のデータ復元部20のハードウエアと同様であってもよいし、共通化されていてもよい。
【0054】
次に、本発明の実施の形態に係る磁気共鳴画像高速再構成法(磁気共鳴イメージング装置1の動作)について説明する。
【0055】
図10に示すように、まず、学習部22が、教師データを用いた学習を行う(ステップS1;深層学習ステップ)。上述のように、学習部22は、高解像度画像25のk空間データの教師データ50(
図6参照)と、教師データ50が間引きされて生成されたk空間データ15を用いて、MRCNN21を学習させる深層学習を行い、MRCNN21の各種パラメータを最適化する。なお、この深層学習ステップでは、所定の空間周波数よりも低周波成分のデータ全体を残し、所定の空間周波数よりも高周波成分についてはランダムにデータを残すアンダーサンプリングフィルタ40を用いて、教師データ50が間引きされ、k空間データ15が生成される。また、この深層学習ステップでは、教師データ50と、MRCNN21を用いて復元された高解像度画像25に対応するk空間データとの誤差に基づいて、誤差逆伝播法による学習が行われる。
【0056】
続いて、データ取得部10は、上記シーケンスに従って動作し、磁気共鳴により観測される核磁気共鳴信号(原信号)に基づくk空間データ15を取得する(ステップS2;データ取得ステップ)。これにより、被検体の断層画像を得るためのk空間データ15が取得される。なお、このデータ取得においては、
図7に示すアンダーサンプリングフィルタ40に従って、シーケンスの繰り返し数又は観測するスライス位置Zの数が規定されている。すなわち、k空間における低周波数成分についてはデータが間引きせずにサンプリングされ、高周波成分についてはランダムに間引きされた状態でデータがサンプリングされる。これにより、データ取得時間(撮像時間)が短縮されている。
【0057】
続いて、データ復元部20が、多重解像度畳み込みニューラルネットワーク(MRCNN21)を用いて、磁気共鳴により観測される原信号に基づくk空間データ15から高解像度画像25を復元する(ステップS3;データ復元ステップ)。具体的には、MRCNN21は、
図4に示すように、k空間データ15から、サイズが異なる(スケールダウンした)k空間データ(1/2データ、1/4データ)を生成する。そして、MRCNN21は、k空間データ15、1/2データ及び1/4データ各々から逆フーリエ変換される画像を、SRCNN35を用いて復元する。さらに、MRCNN21は、復元されたデータを、アップサンプリング、すなわちデータのサイズを合わせて重ね合わせることにより、高解像度画像25を復元する。
【0058】
図11には、劣化した被検体の断層画像が示され、
図12には、本実施の形態に係る磁気共鳴画像高速再構成法において復元された被検体の断層画像が示されている。
図11及び
図12を比較するとわかるように、被検体の断層画像が高精細に復元されている。なお、この断層画像は、T1強調画像である。T1強調画像は、2次元で脳の組織を見るための画像である。T1強調画像は、大脳皮質又は骨などの解剖学的な構造を見る時に使用される。T1は、縦磁化の緩和時間である。T1強調画像は、RFパルスが出射される間隔である繰り返し時間を短くしてT1の影響を強調し、RFパルスを出射してからエコー信号の中央までの時間であるエコー時間を短くしてT2(横磁化の緩和時間)の影響を少なくすることにより得られる画像である。
【0059】
本実施の形態に係る磁気共鳴画像高速再構成法による画像の復元状態を評価した。データとしては、頭部のMRI断層画像(サイズ:256×256)を用いた。trainを100枚、testを50枚として、本実施の形態に係るMRCNN21によるデータ復元、SRCNN35単独でのデータ復元との間での性能比較を行った。画像の復元の評価方法としては、PSNR、MISSIM、KPSNRが用いられた。
【0060】
まず、PSNRについて説明する。xを原画像、yを比較したい画像とすると、PSNRは、以下の式で定義される。
【数15】
【数16】
【0061】
上式により、PSNRの値が大きいほど原画像に近いという評価になる。
SSIMは、画像の小領域毎に計算される。計算式を以下に示す。
【数17】
ただし、
【数18】
ここで、x,yは、それぞれの小領域の画素であり、μ
x,μ
yは、それぞれの平均画素値であり、σx,σyは、それぞれの標準偏差であり、σ
xyはxとyの共分散である。また、C
1=(K
1L)
2
,C
2=(K
2L)
2であり、Lはダイナミックレンジ、K1、K2は、任意の値であり、例えば、K
1=0.01、K
2=0.03とすることができる。
【0062】
ここで、α=β=γ=1かつC
3=C
2/2とすると、SSIM(x,y)は以下の式で表される。
【数19】
SSIMは小領域ごとに計算され、画像に対するSSIMの平均を取ったものをMSSIMという。MSSIMを次式に表す。
【数20】
【0063】
次に、KPSNRについて説明する。KPSNRは周波数領域ごとの復元度を調べるために考案した手法で、K空間でPSNRを測定する。
図13、
図14、
図15の白の領域をそれぞれ低周波領域、中間領域、高周波領域として抽出し、それぞれの領域でKPSNRを測定した。各領域での測定結果を、KPSNR(low)、KPSNR(middle)、KPSNR(high)とした。
【0064】
図16に示すように、MRCNN21は、PSNR、MSSIM、KPSNR(low)、KPSNR(middle)、KPSNR(high)全てにおいて、SRCNN35単独の場合よりも、評価値が良好であった。すなわち、畳み込みニューラルネットワークの多重化を行った効果、本実施の形態に係るMRCNN21の効果を確認することができた。
【0065】
また、本実施の形態に係る磁気共鳴画像高速再構成法を、MRA(Magnetic Resonance Angiography)画像に適用した。MRA画像とは、3次元の血管の構造を把握するために使用される画像である。MRA画像は、MIP(Maximum Intensity Projection;最大投影法)処理を行うことにより生成され、血管が強調された画像となる。
【0066】
MIP処理とは、3次元の画像データを任意の方向(投影方向)にある2次元の投影面に投影した場合に、3次元の画像データにおける投影面の同じ位置に投影される複数の画素の最大値を投影面上の各画素の値とする処理のことである。この処理により、3次元の画像データにおいて投影方向に沿った複数のボクセルの中で最も信号強度が高いボクセルの信号強度が、MRA画像のピクセルの信号強度となる。
【0067】
ここでは、患者1人当たりのMRA画像のデータのサイズを、175×512×512(スライス×縦×横)とした。また、訓練データを200枚とし、テストデータを100枚とした。フィルタとしては、
図7のアンダーサンプリングフィルタを用いた。また、学習時のハイパーパラメータとして、学習係数(誤差逆伝播の係数)を1.0×10
-5とし、反復回数を1.0×10
4とし、バッチサイズを10とした。
【0068】
図17(A)~
図17(C)には、頭頂から見た場合の、すなわち頭頂から見た場合に視線方向に投影されるMIP画像の一例が示されている。
図17(A)は、元画像(正解画像)である。この正解画像は、時間をかけて計測された高精細な画像である。ここで、この正解画像を教師データ50としてMRCNN21により学習を行った。例えば、アンダーサンプリングフィルタ40を用いて
図17(A)に示す元画像から、
図17(C)に示す劣化画像を生成し、MRCNN21の学習を行った。さらに、学習されたMRCNN21により、
図17(C)に示す劣化画像から、
図17(B)に示す復元画像が得られた。
【0069】
図18(A)~(C)には、
図17(A)~(C)の白枠部分の拡大図が示されている。
図18(A)と
図18(B)とを比較して示すように、MRCNN21で復元された復元画像は、元画像(正解画像)と比べ、精細さに遜色はなかった。
【0070】
図19(A)~
図19(C)には、横方向から見た場合のMIP画像の一例が示されている。
図19(A)は、元画像(正解画像)である。この正解画像は、時間をかけて高精細な画像である。同様に、この正解画像を教師データ50としてMRCNN21により学習を行った。例えば、アンダーサンプリングフィルタ40を用いて
図17(A)に示す元画像から、
図19(C)に示す劣化画像を生成し、MRCNN21の学習を行った。そして、学習されたMRCNN21により、
図19(C)に示す劣化画像から、
図19(B)に示す復元画像が得られた。
【0071】
図20(A)~(C)には、
図19(A)~(C)の白枠部分の拡大図が示されている。
図20(A)と
図20(B)とを比較して示すように、MRCNN21で復元された復元画像は、元画像(正解画像)と比べ、精細さに遜色はなかった。
【0072】
図21(A)~
図21(C)には、正面方向から見た場合のMIP画像の一例が示されている。
図21(A)は、元画像(正解画像)である。この正解画像は、時間をかけて高精細な画像である。ここで、この正解画像を教師データ50としてMRCNN21により学習を行った。例えば、アンダーサンプリングフィルタ40を用いて
図21(A)に示す元画像から、
図21(C)に示す劣化画像を生成し、MRCNN21の学習を行った。そして、学習されたMRCNN21により、
図21(C)に示す劣化画像から、
図21(B)に示す復元画像が得られた。
【0073】
図22(A)~(C)には、
図21(A)~(C)の白枠部分の拡大図が示されている。
図22(A)と
図22(B)とを比較して示すように、MRCNN21で復元された復元画像は、元画像(正解画像)と比べ、精細さに遜色はなかった。
【0074】
なお、画像の復元度の評価を行うために、上記式(14)のPSNRが使用された。
図23には、PSNRの結果が示されている。具体的には、MIP処理を行う前のMRA原原画像について、元画像に対する劣化画像のPSNR値と、元画像に対する復元画像のPSNR値と、復元画像のPSNR値とが求められ、それらの差が復元度として求められた。MIP画像(頭頂)と、MIP画像(横方向)と、MIP画像(正面方向)とについても、同様にPSNR値が求められた。
図23に示すように、MIP画像(頭頂)と、MIP画像(横方向)と、MIP画像(正面方向)とについて、それぞれの復元度は大幅に向上しており、MIP処理前のMRA原画像と、MIP処理後の画像とを比較すると、MIP処理を行った方がより復元度が向上していることが分かった。これは、MIP処理後の画像は、投影方向に沿ったボクセルの最大値に基づく画像であり、MIP処理前に比べノイズ成分が低減されているためであると考えられる。
【0075】
以上詳細に説明したように、本実施の形態によれば、処理時間が圧縮センシングよりも短い多重解像度畳み込みニューラルネットワーク(MRCNN21)を用いて観測点数の少ない第1のk空間データ15をスケールダウンし、多重化して超解像処理を行って得られる画像を用いて画像の再構成を行う。このため、観測点数を少なくして撮像時間を短縮しつつ、多重解像度畳込みニューラルネットワークを用いて再構成時間を短縮して、高周波成分だけでなく低周波成分も損なうことなく高解像度画像25を復元することができる。この結果、撮像時間及び再構成時間を双方短縮しつつ、高解像度画像25を復元することができる。
【0076】
MRI断層撮影では、周波数領域でアンダーサンプリングが行われるため、様々な周波数成分が不足してしまう。また、SRCNN35単体では、フィルタのサイズが小さいため、高周波成分しかとらえることができないという問題があり、フィルタを大きくするしかなく、この場合にはパラメータ数が増えて学習が困難になる。そこで、本実施の形態に係る磁気共鳴イメージング装置1では、超解像畳み込みニューラルネットワークの多重化を行って、入力画像のサイズを小さくすることで相対的にフィルタのサイズを大きくしたので、低周波成分の復元も行うことができるようになり、この結果、高精細な画像を復元することができるようになった。
【0077】
また、本実施の形態では、画像ではなく、磁気共鳴信号(原信号)から構成されるk空間そのものに対する学習により、データ復元を行う。これにより、低域から高域までの様々な空間周波数成分を直接復元することができるようになるので、より高精細な画像の復元が可能となる。このため、これまで発見が困難であった病変の発見も可能になることが期待される。
【0078】
また、本実施の形態では、画像の復元に有益な低周波成分全体を残し、病変を発見する上で重要な高周波成分をランダムに残すアンダーサンプリングフィルタを用いて学習を行った。これにより、病変を発見し易く、かつ、高精細な画像を復元することが可能となった。
【0079】
なお、本実施の形態では、アンダーサンプリングフィルタ40のデータ取得率を10%としたが、本発明はこれには限られない。データ取得率は、10%を超えてもよいし、10%未満であってもよい。
【0080】
なお、アンダーサンプリングフィルタ40の形状は、
図7に示すものには限られない。アンダーサンプリングフィルタ40の形状は、原信号の低周波成分と、高周波成分とが大きく損なわれないようなものであるのが望ましい。
【0081】
また、学習部22は、高解像度k空間データ15と教師データ50との誤差に基づく誤差逆伝播法により学習を行った。これにより、k空間データにおける誤差に基づいて直接学習が行われるようになるので、学習時間を短縮することが可能となる。なお、学習方法は、誤差逆伝播法には限られない。MRCNN21のパラメータを最適化できるのであれば、様々な学習方法を適用可能である。
【0082】
なお、上記実施の形態では、MRCNN21における多重化を3つとしたが、2つでもよいし、4つ以上でもよい。多重化が2つの場合には、1/2データが生成され、多重化が4つ以上の場合には、1/2データ、1/4データ、1/8データ、・・・を生成して、それぞれをSRCNN35に入力すればよい。
【0083】
また、多重化するのはSRCNN35に限られない。他の畳み込みニューラルネットワークを多重化するようにしてもよい。
【0084】
なお、本発明は、データ取得(撮影)のパルスシーケンスに制限はない。スピンエコー法、高速スピンエコー法、グラディエントエコー法、エコープレナー法、ハーンエコー法、スティミュレイテドエコー法など、様々な方法に対応することができる。
【0085】
その他、データ復元部20のハードウエア構成やソフトウエア構成は一例であり、任意に変更および修正が可能である。
【0086】
制御部61、主記憶部62、外部記憶部63、操作部64、表示部65及び入出力部66、内部バス60などから構成されるデータ復元部20の処理を行う中心となる部分は、上述のように、専用のシステムによらず、通常のコンピュータシステムを用いて実現可能である。例えば、前記の動作を実行するためのコンピュータプログラムを、コンピュータが読み取り可能な記録媒体(フレキシブルディスク、CD-ROM、DVD-ROM等)に格納して配布し、当該コンピュータプログラムをコンピュータにインストールすることにより、前記の処理を実行するデータ復元部20を構成してもよい。また、インターネット等の通信ネットワーク上のサーバ装置が有する記憶装置に当該コンピュータプログラムを格納しておき、通常のコンピュータシステムがダウンロード等することで磁気共鳴イメージング装置1を構成してもよい。
【0087】
コンピュータの機能を、OS(オペレーティングシステム)とアプリケーションプログラムの分担、またはOSとアプリケーションプログラムとの協働により実現する場合などには、アプリケーションプログラム部分のみを記録媒体や記憶装置に格納してもよい。
【0088】
搬送波にコンピュータプログラムを重畳し、通信ネットワークを介して配信することも可能である。たとえば、通信ネットワーク上の掲示板(BBS, Bulletin Board System)にコンピュータプログラムを掲示し、ネットワークを介してコンピュータプログラムを配信してもよい。そして、このコンピュータプログラムを起動し、OSの制御下で、他のアプリケーションプログラムと同様に実行することにより、前記の処理を実行できるように構成してもよい。
【0089】
この発明は、この発明の広義の精神と範囲を逸脱することなく、様々な実施の形態及び変形が可能とされるものである。また、上述した実施の形態は、この発明を説明するためのものであり、この発明の範囲を限定するものではない。すなわち、この発明の範囲は、実施の形態ではなく、特許請求の範囲によって示される。そして、特許請求の範囲内及びそれと同等の発明の意義の範囲内で施される様々な変形が、この発明の範囲内とみなされる。
【産業上の利用可能性】
【0090】
本発明は、磁気共鳴画像の撮影時間、データ処理時間の短縮及び高解像度化に適用することができる。
【符号の説明】
【0091】
1 磁気共鳴イメージング装置、10 データ取得部、15 k空間データ、20 データ復元部、21 多重解像度畳み込みニューラルネットワーク(MRCNN)、22 学習部、25 高解像度画像、35 超解像畳み込みニューラルネットワーク(SRCNN)、40 アンダーサンプリングフィルタ、41 フーリエ変換部、45 パラメータ設定部、50 教師データ、60 内部バス、61 制御部、62 主記憶部、63 外部記憶部、64 操作部、65 表示部、66 入出力部、69 プログラム