JP2020146449A

JP2020146449A - 磁気共鳴画像高速再構成法及び磁気共鳴イメージング装置

Info

Publication number: JP2020146449A
Application number: JP2020025707A
Authority: JP
Inventors: 竹内　純一; Junichi Takeuchi; 純一竹内; 豊實松; Yutaka JITSUMATSU; 雅則川喜田; Masanori Kawakita; 自然北崎; Shizen Kitazaki; 久原　重英; Shigehide Kuhara; 重英久原
Original assignee: Kyushu University NUC; Kyorin University
Current assignee: Kyushu University NUC; Kyorin University
Priority date: 2019-03-06
Filing date: 2020-02-18
Publication date: 2020-09-17
Anticipated expiration: 2040-02-18
Also published as: JP7466141B2

Abstract

【課題】撮像時間及び再構成時間を双方短縮しつつ、高解像度画像を復元することができる磁気共鳴画像高速再構成法及び磁気共鳴イメージング装置を提供する。【解決手段】磁気共鳴画像高速再構成法は、多重解像度畳み込みニューラルネットワークを用いて、磁気共鳴により観測される原信号に基づく第１のｋ空間データから高解像度画像を復元するデータ復元ステップ（ステップＳ３）を含む。多重解像度畳み込みニューラルネットワークでは、第１のｋ空間データから、サイズが異なる第２のｋ空間データを生成し、第１のｋ空間データ及び第２のｋ空間データ各々から得られる画像を、超解像畳み込みニューラルネットワークを用いて超解像化し、超解像化された画像を、サイズを合わせて重ね合わせることにより、高解像度画像を復元する。【選択図】図１０

Description

本発明は、磁気共鳴画像高速再構成法及び磁気共鳴イメージング装置に関する。

磁気共鳴画像法（Magnetic Resonance Imaging; MRI）は生体内の情報を得る有力な手段の１つである。しかし、ＭＲＩ装置の検査には、数十分から１時間ほどの時間が必要となる。そこで、患者の負担軽減のため、撮像時間のさらなる短縮が求められている。この１０年間、情報処理によるＭＲＩの高速化と高精細化に貢献してきたのが圧縮センシング（Compressed Sensing; CS）である（例えば、特許文献１、非特許文献１参照）。圧縮センシングは、アンダーサンプリングを行って撮像時間を削減しつつ、計算により少ない観測データから元の画像を復元し再構成する技術である。

特開２０１６−１２３８５３号公報

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.

しかし、圧縮センシングのアルゴリズムは観測数の３乗オーダーの計算量を必要とする。このため、観測データを少なくして削減した撮像時間よりも、画像を復元し再構成する再構成時間の方が超えてしまうといった不都合が生ずる場合がある。

本発明は、上記実情の下になされたものであり、撮像時間及び再構成時間を双方短縮しつつ、高解像度画像を復元することができる磁気共鳴画像高速再構成法及び磁気共鳴イメージング装置を提供することを目的とする。

上記目的を達成するために、本発明の第１の観点に係る磁気共鳴画像高速再構成法は、
多重解像度畳み込みニューラルネットワークを用いて、磁気共鳴により観測される原信号に基づく第１のｋ空間データから高解像度画像を復元するデータ復元ステップを含み、
前記多重解像度畳み込みニューラルネットワークでは、
前記第１のｋ空間データから、サイズが異なる第２のｋ空間データを生成し、
前記第１のｋ空間データ及び前記第２のｋ空間データ各々から得られる画像を、超解像畳み込みニューラルネットワークを用いて超解像化し、
超解像化された画像を、サイズを合わせて重ね合わせることにより、前記高解像度画像を復元する。

この場合、前記データ復元ステップに先立って、
高精細なｋ空間データを教師データとし、前記教師データが間引きされて生成されたｋ空間データを前記第１のｋ空間データとして、前記多重解像度畳み込みニューラルネットワークを学習させる深層学習ステップを含む、
こととしてもよい。

また、前記深層学習ステップでは、
所定の空間周波数よりも低周波成分のデータ全体を残し、前記空間周波数よりも高周波成分についてはランダムにデータを残すアンダーサンプリングフィルタを用いて、前記教師データを間引きすることにより、学習用の前記第１のｋ空間データを生成する、
こととしてもよい。

前記深層学習ステップでは、
前記教師データと、前記多重解像度畳み込みニューラルネットワークを用いて復元された前記高解像度画像に対応する第３のｋ空間データとの誤差に基づく誤差逆伝播法を用いて学習を行う、
こととしてもよい。

本発明の第２の観点に係る磁気共鳴イメージング装置は、
磁気共鳴により観測される原信号に基づく第１のｋ空間データを取得するデータ取得部と、
多重解像度畳み込みニューラルネットワークを用いて、前記第１のｋ空間データから高解像度画像を復元するデータ復元部と、
を備え、
前記多重解像度畳み込みニューラルネットワークは、
前記第１のｋ空間データから、サイズが異なる第２のｋ空間データを生成し、
前記第１のｋ空間データ及び前記第２のｋ空間データ各々から得られる画像を、超解像畳み込みニューラルネットワークを用いて超解像化し、
超解像化された画像を、サイズを合わせて重ね合わせることにより、前記高解像度画像を復元するように構成されている。

本発明によれば、処理時間が圧縮センシングよりも短い多重解像度畳み込みニューラルネットワークを用いて観測点数の少ない第１のｋ空間データをスケールダウンし、多重化して超解像処理を行って得られる画像を用いて画像の再構成を行う。このため、観測点数を少なくして撮像時間を短縮しつつ、多重解像度畳込みニューラルネットワークを用いて再構成時間を短縮して、高周波成分だけでなく低周波成分も損なうことなく画像を復元することができる。この結果、撮像時間及び再構成時間を双方短縮しつつ、高解像度画像を復元することができる。

本発明の実施の形態に係る磁気共鳴イメージング装置の構成を示すブロック図である。（Ａ）は、被検体に印加されるＸＹ方向の傾斜磁場を示す図である。（Ｂ）は、Ｚ軸方向の傾斜磁場を示す図である。（Ａ）は、２次元のｋ空間データを示す図である。（Ｂ）は、スライス画像を示す図である。多重解像度畳み込みニューラルネットワークの構成を示すブロック図である。超解像畳み込みニューラルネットワークの構成を示す模式図である。学習部の構成を示すブロック図である。アンダーサンプリングフィルタの一例を示す図である。誤差逆伝搬法を用いた深層学習によるパラメータの更新の様子を示す模式図である。データ復元部のハードウエアの構成を示すブロック図である。本発明の実施の形態に係る磁気共鳴画像高速再構成法を示すフローチャートである。データ取得部で取得されたｋ空間データから再構成された断層画像（Ｔ１強調画像）を示す図である。データ復元部で復元されたｋ空間データから再構成された断層画像（Ｔ１強調画像）を示す図である。ＰＳＮＲを測定する低周波領域を示す模式図である。ＰＳＮＲを測定する中間領域を示す模式図である。ＰＳＮＲを測定する高周波領域を示す模式図である。Ｔ１強調画像の復元の評価結果を示す図である。（Ａ）は、頭頂方向から見たときに得られるＭＩＰ画像の元画像の一例である。（Ｂ）は、頭頂方向から見たときに得られるＭＩＰ画像の多重解像度畳み込みニューラルネットワークによる復元画像である。（Ｃ）は、頭頂方向から見たときに得られるＭＩＰ画像をアンダーサンプリングフィルタにより劣化させた劣化画像である。（Ａ）は、図１７（Ａ）の画像の白枠部分である。（Ｂ）は、図１７（Ｂ）の画像の白枠部分である。（Ｃ）は、図１７（Ｃ）の画像の白枠部分である。（Ａ）は、横方向から見たときに得られるＭＩＰ画像の元画像の一例である。（Ｂ）は、横方向から見たときに得られるＭＩＰ画像の多重解像度畳み込みニューラルネットワークによる復元画像である。（Ｃ）は、横方向から見たときに得られるＭＩＰ画像をアンダーサンプリングフィルタにより劣化させた劣化画像である。（Ａ）は、図１９（Ａ）の画像の白枠部分である。（Ｂ）は、図１９（Ｂ）の画像の白枠部分である。（Ｃ）は、図１９（Ｃ）の画像の白枠部分である。（Ａ）は、正面方向から見たときに得られるＭＩＰ画像の元画像の一例である。（Ｂ）は、正面方向から見たときに得られるＭＩＰ画像の多重解像度畳み込みニューラルネットワークによる復元画像である。（Ｃ）は、正面方向から見たときに得られるＭＩＰ画像をアンダーサンプリングフィルタにより劣化させた劣化画像である。（Ａ）は、図２１（Ａ）の画像の白枠部分である。（Ｂ）は、図２１（Ｂ）の画像の白枠部分である。（Ｃ）は、図２１（Ｃ）の画像の白枠部分である。劣化画像から復元画像への復元度を比較する表である。

以下、本発明の実施の形態について図面を参照して詳細に説明する。各図面においては、同一又は同等の部分に同一の符号を付す。

図１に示すように、磁気共鳴イメージング装置１は、データ取得部１０と、データ復元部２０と、を備える。

データ取得部１０は、静磁場の中に置かれた被検体（生体）に一連のシーケンスに従って高周波磁場及び傾斜磁場を印加して、被検体から発生する核磁気共鳴信号（原信号）を検出する。核磁気共鳴信号は、被検体内の水素密度分布を示す空間周波数の領域（ｋ空間）の信号として観測される。データ取得部１０は、磁気共鳴により観測される原信号に基づく第１のｋ空間データとしてのｋ空間データ１５を取得する。

具体的には、データ取得部１０は、図２（Ａ）に示すように、被検体を＋Ｚ方向の静磁場Ｂ_０に置いた状態で、静磁場に垂直な方向から９０°パルス及び１８０°パルスなどのＲＦ（Radio Frequency）パルスを被検体に印加する。これにより、被検体の水素原子核が励起されて、核磁気共鳴信号（エコー信号）を取り出すことができる。データ取得部１０は、この核磁気共鳴信号（エコー信号）をサンプリングして取得する。データ取得部１０は、ｋ空間データ１５を、サンプリング点を間引きつつ取得することができる。

データ取得部１０において、ＲＦパルスは、図２（Ｂ）に示すように、Ｚ軸方向に勾配を有するスライス傾斜磁場Ｂ（Ｚ）と同時に被検体に印加される。通常は、Ｚ軸方向が被検体の体軸方向となる。スライス位置Ｚに対応する静磁場強度をＢ（Ｚ）とすると、スライス位置Ｚ_０を選択励起するため、周波数γＢ（Ｚ_０）を持つ正弦波のＲＦパルスを印加する。ここで、γは磁気回転比である。

さらに、図２（Ａ）に示すように、データ取得部１０において、被検体には、Ｙ軸方向に傾きＧｙを有する位相エンコード傾斜磁場と、Ｘ軸方向に傾きＧｘを有する周波数エンコード傾斜磁場とが加えられる。さらに、核磁気共鳴信号（エコー信号）を観測するタイミングで、Ｘ軸方向に傾きＧｘを有する周波数エンコード傾斜磁場が再度印加される。上述のようなシーケンスが位相エンコード傾斜磁場の強度（傾きＧｙ）を変えながら繰り返され、その度に核磁気共鳴信号（エコー信号）が観測される。

上記シーケンスにおいて、観測される複数の核磁気共鳴信号（エコー信号）により、２次元のｋ空間データが生成される。図３（Ａ）に示すように、観測される核磁気共鳴信号（エコー信号）は、ｋ空間では、点線で示される空間周波数Ｋｙが同じ直線の領域に対応するデータとなる。

ｋ空間における位置座標（Ｋｘ，Ｋｙ）は、周波数エンコード傾斜磁場の傾きＧｘ、位相エンコード傾斜磁場の傾きＧｙで表現すれば（γＧｘｔ，γＧｙｔ）となる。ここで、ｔは、位相エンコード傾斜磁場の傾きＧｙ、周波数エンコード傾斜磁場の傾きＧｘの印加時間である。核磁気共鳴信号（エコー信号）は、位相エンコード傾斜磁場の傾きＧｙの強度を変えながら観測されるため、それぞれの信号はＫｙ軸の位置が異なるデータとなる。このようにして、スライス位置Ｚ_０における２次元のｋ空間データが構成される。

データ取得部１０は、このような計測を、ＲＦパルスの周波数を変更しながら、すなわちスライス位置Ｚを変更しながら繰り返し行う。これにより、スライス位置Ｚが異なる２次元ｋ空間データが観測される。なお、このｋ空間データを逆フーリエ変換すれば、そのＺ位置におけるＸＹ面の２次元画像が得られる。その画像データを図３（Ｂ）に示すように並べれば、被検体の３次元の画像データが得られるようになる。

なお、データ取得部１０でデータが取得される時間、すなわち撮像時間を短縮する場合には、シーケンスの繰り返し数又は観測するスライス位置Ｚの数を少なくする必要がある。シーケンスの繰り返し数を少なくすれば、図３（Ａ）に点線で示されるエコー信号の本数が少なくなり、スライス位置Ｚの数を少なくすれば、図３（Ｂ）に示される２次元画像の枚数が少なくなる。

図１に戻り、データ復元部２０は、データ取得部１０で取得されたｋ空間データ１５を入力する。データ復元部２０は、多重解像度畳み込みニューラルネットワークを用いて、ｋ空間データ１５から高解像度画像２５を復元する。また、畳み込みニューラルネットワークは、もともと人間の脳の視覚に関わる構造を参考にして考案され、名前の通り畳み込み層を持つことを特徴とする。

データ復元部２０は、多重解像度畳み込みニューラルネットワーク２１を用いてｋ空間データ１５を高解像化し、高解像度画像２５を生成して出力する。多重解像度畳み込みニューラルネットワーク２１は、多重解像度（Multi-Resolution）を、畳み込みニューラルネットワーク（Convolutional Neural Network; ＣＮＮ）で実現していることから、ＭＲＣＮＮ２１ともいう。

ＭＲＣＮＮ２１は、入力したｋ空間データ１５から、サイズの異なる（スケールダウンした）ｋ空間データ（第２のｋ空間データ）を複数生成する。ＭＲＣＮＮ２１は、図４に示すように、まず、ｋ空間データ１５を１／２、１／４にリサイズして低周波成分を取得し、２つのｋ空間データ（１／２データ、１／４データ）を生成する。ｋ空間データ１５、１／２データ及び１／４データはそれぞれ逆フーリエ変換される。例えば、入力されたｋ空間データ１５が逆フーリエ変換されて生成される入力画像が、２５６×２５６ビットのデータであった場合には、１２８×１２８ビットの画像と、６４×６４ビットの画像とが生成される。このリサイズ処理には、例えばバイキュービック補間が用いられる。

続いて、サイズの異なる画像データそれぞれに対して超解像畳み込みニューラルネットワーク（ＳＲＣＮＮ）３５を用いた画像復元を行う。超解像とは、与えられた低解像度画像からより解像度の高い高解像度画像を得る手法のことである。ここでは、図４に示すように、入力画像、１／２データに対応する画像、１／４データに対応する画像に対してそれぞれＳＲＣＮＮ３５により超解像処理が行われる。

ＳＲＣＮＮ３５は、図５に示す３層の畳み込み層で構成され、大きく分けて３つのステップを含んでいる。第１のステップは、パッチ分解・抽出であり、第２のステップは、辞書学習・パッチ推定（非線形写像）であり、第３のステップは、データの再構成である。

（第１のステップ）
第１のステップでは、パッチ分解・特徴量の抽出などが行われる。第１のステップでの写像（特徴マップ）を以下の関数Ｆ_１で表す。パッチとは画像を細かく分解したものである。パッチ分解することで２次元画像データのサイズが小さくなる。
ここで、Ｗ_１は畳み込みのためのフィルタ群（重みフィルタ）であり、Ｂ_１は、バイアスである。Ｗ_１のサイズをｃ×ｆ_１×ｆ_１×ｎ_１とする。ｃは入力データのチャネル数であり、ｆ_１は、フィルタの空間的なサイズ（縦×横）を表し、例えば９である。ｎ_１は、特徴マップの数であり、例えば６４である。

（第２のステップ）
第２のステップでは、低解像度パッチを高い解像度パッチに非線形的に写像し、圧縮センシングの計算を行う。データのサイズが小さいため、このときの圧縮センシングについては、短時間のうちに行うことができる。第２のステップである非線形写像に相当する写像（特徴マップ）を関数Ｆ_２で表す。
ここで、Ｗ_２は重みフィルタを表す。Ｗ_２のサイズは、ｎ_１×１×１×ｎ_２である。ｎ_２は、例えば３２である。Ｂ_２は、バイアスを表すｎ_２次元ベクトルである。

第３のステップでは、高解像度パッチベクトルを結合し、推定高解像度データｄとして出力する。第３のステップである写像は、以下の関数Ｆで表すことができる。
ここで、Ｗ_３は重みフィルタを表す。Ｗ_３のサイズは、ｎ_２×ｆ_３×ｆ_３×ｃである。ｆ_３は、例えば５である。Ｂ_３は、バイアスを表すｃ次元ベクトルである。ここで出力されるデータが推定高解像度データｄとなる。

ここで、図４に示すように、ｋ空間データ１５に対応する入力画像が入力されるＳＲＣＮＮ３５で推定される推定高解像度データｄをｄ１とする。また、１／２データに対応する画像が入力されるＳＲＣＮＮ３５で推定される推定高解像度データｄをｄ２とする。さらに、１／４データに対応する画像が入力されるＳＲＣＮＮ３５で推定される推定高解像度データｄをｄ３とする。ＳＲＣＮＮ３５で復元された推定高解像度データｄ３は、２つの畳み込み層（ＣＯＮＶ）を経て、アップサンプリングされ、推定高解像度データｄ２と重ね合わせられる。ここで、アップサンプリングとは、推定高解像度データを拡大する処理であり、アンプーリング（プーリングの逆の処理）とも呼ばれる。例えば、重ね合わせられた推定高解像度データｄ２＋ｄ３は、２つの畳み込み層（ＣＯＮＶ）を経てアップサンプリングされ、ｋ空間データ１５に対応する画像から生成された推定高解像度データｄ１と重ね合わせられる。この重ね合わされた推定高解像度データｄ１＋ｄ２＋ｄ３が、２つの畳み込み層（ＣＯＮＶ）を経て、復元された高解像度画像として復元される。この高解像度画像２５は、学習時にフーリエ変換され、ｋ空間データとして出力される。なお、各段の２つの畳み込み層のフィルタサイズは、３×３×２×２又は３×３×２×１である。

データ復元部２０において、ＭＲＣＮＮ２１によりｋ空間データ１５から得られる画像の高解像化を行うためには、ＭＲＣＮＮ２１のパラメータを最適化する必要がある。このため、データ復元部２０は、深層学習により、ＭＲＣＮＮ２１を学習させ、そのパラメータを最適化する学習部２２を備えている。

学習部２２は、教師データ５０を有している。教師データ５０としては、間引きすることなくサンプリングされた被検体の高精細な断層画像をフーリエ変換して得られるｋ空間データが用いられる。学習部２２は、教師データ５０の一部が間引きされて生成されたデータを、ｋ空間データ１５としてＭＲＣＮＮ２１に入力する。

学習部２２は、アンダーサンプリングフィルタ４０を有している。アンダーサンプリングフィルタ４０は、データからサンプルを間引いて解像度を下げるフィルタである。図７では、ｋ空間において、白い部分がデータを取得する部分であり、黒い部分がデータを間引く部分である。図７に示すように、アンダーサンプリングフィルタ４０は、ｋ空間データにおける中間の空間周波数よりも低周波成分（白い領域）のデータ全体を残しつつ、空間周波数よりも高周波成分についてはランダムにデータを残すような形状のフィルタとなっている。アンダーサンプリングフィルタ４０は、全体のデータ取得率を１０％としている。例えば、低周波成分と高周波成分とのデータ取得率を１：１とすることができるが、この割合に特に限定されるものではない。

このアンダーサンプリングフィルタ４０は、低周波成分にエネルギーのほとんどが中心に集中しており、画像全体を高精度に復元することができるように構成されている。また、アンダーサンプリングフィルタ４０は、また、高周波成分は、病変を発見する上で重要であるため、ポアソンディスクサンプリングによりランダムに高周波成分を取得するように構成されている。

学習部２２は、アンダーサンプリングフィルタ４０を通過して間引きされたｋ空間データ１５をＭＲＣＮＮ２１に入力する。ＭＲＣＮＮ２１は、高解像度画像２５を出力する。学習部２２は、この場合に、ＭＲＣＮＮ２１で生成される高解像度画像２５がフーリエ変換部４１でフーリエ変換されて生成されるｋ空間データが教師データ５０に近づくようにＭＲＣＮＮ２１のパラメータ調整を行う。

学習部２２は、パラメータ設定部４５を備えている。パラメータ設定部４５は、教師データ５０と、教師データ５０が間引きされて生成されたｋ空間データ１５から復元された高解像度画像２５に対応するｋ空間データ（第３のｋ空間データ）との誤差に基づいて、誤差逆伝搬法を用いて深層学習を行う。誤差逆伝搬法は、入力層、中間層、出力層からなるネットワークに対し、出力層から入力層にかけて誤差の勾配を逆伝搬させることで各層の重みフィルタＷ_１，Ｗ_２，Ｗ_３とバイアスＢ_１，Ｂ_２，Ｂ_３といったパラメータを更新する教師付き学習アルゴリズムである。

各層の重みについてパラメータの更新の様子を説明する。出力層付近に着目したネットワークの構造を図８に示す。Ｎ層第ｋニューロンへの入力、Ｎ層第ｋニューロンの出力をそれぞれｘ_ｋ ^（Ｎ），ｙ_ｋ ^（Ｎ）とする。また、Ｎ−１層第ｊニューロンからＮ層第ｋニューロンへの重みをｗ_ｊｋ ^（Ｎ）とする。Ｎ層第ｋニューロンの活性化関数をｇ_ｋ ^（Ｎ）とおくと、ｘ_ｋ ^（Ｎ），ｙ_ｋ ^（Ｎ）は、ニューロンモデルから次の式で表すことができる。

あるネットワークの入力に対する出力層の出力（高解像度画像２５から得られるｋ空間データ）と教師データ５０との誤差を表す関数を誤差関数Ｌ（ｙ_１ ^（Ｎ），・・・，ｙ_ｋ ^（Ｎ），・・・，ｙ_Ｋ ^（Ｎ））とする。まず、誤差関数ＬのＮ−１層第ｊニューロンからＮ層第ｋニューロンへの重みｗ_ｊｋ ^（Ｎ）に関する勾配は、式（４）、（５）と連鎖律によって式（６）で表される。
式（４）より、以下の式が成立する。

ただし、
のように定義している。このように、誤差関数Ｌの勾配は、式（７）に示すように、出力側の層による項（式（８））と、出力層の一層前の出力ｙ_ｊ ^{（Ｎ−１）}から計算することができる。

次に、Ｎ−１層（中間層）の重みに関する勾配を考える。同様に、誤差関数ＬのＮ−２層第ｉニューロンからＮ−１層第ｊニューロンへの重みｗ_ｉｊ ^{（Ｎ−１）}に関する勾配は、次式で表される。
ここで、
は、出力層から逆伝播された
から計算することができる。このように誤差関数を入力層に向かって逆伝播させていき、各層のパラメータに関する誤差関数の勾配を求めていく方法が、誤差逆伝搬法である。

パラメータの更新は、最急降下法の考え方に基づいて行う。最急降下法では誤差関数Ｌを用いて次式のように各パラメータを更新する。ただし、ｔは、反復回数を表す。
重みの更新量は、誤差関数Ｌの勾配に学習係数ηをかけたものを重みの更新量とする。

最急降下法では学習率が低いと収束が遅くなってしまう。また、ｔに対して勾配が振動するような場合も収束が遅くなる。これを解決するためにモーメンタム法と呼ばれる手法を用いるようにしてもよい。モーメンタム法における更新式を以下に示す。
μをモーメンタムと呼ぶ。上式（１２）は、一回前の反復時の情報を用いて勾配の振動を抑制している。また勾配が振動せず同じ方向の場合は更新量を増加させる役割を持つ。なお、勾配法としてはａｄａｍ法を用いるようにしてもよい。

図９に示すように、データ復元部２０は、制御部６１、主記憶部６２、外部記憶部６３、操作部６４をハードウエア構成として備えている。主記憶部６２、外部記憶部６３、操作部６４はいずれも内部バス６０を介して制御部６１に接続されている。

制御部６１は、ＣＰＵ（Central Processing Unit）等から構成されている。このＣＰＵが、外部記憶部６３に記憶されているプログラム６９に従ってＭＲＣＮＮ２１及び学習部２２の処理を実行することにより、図１に示すデータ復元部２０が実現される。

主記憶部６２は、ＲＡＭ（Random-Access Memory）等から構成されている。主記憶部６２には、外部記憶部６３に記憶されているプログラム６９がロードされる。この他、主記憶部６２は、制御部６１の作業領域（データの一時記憶領域）として用いられる。

外部記憶部６３は、フラッシュメモリ、ハードディスク、ＤＶＤ−ＲＡＭ（Digital Versatile Disc Random-Access Memory）、ＤＶＤ−ＲＷ（Digital Versatile Disc ReWritable）等の不揮発性メモリから構成される。外部記憶部６３には、制御部６１に実行させるためのプログラム６９があらかじめ記憶されている。また、外部記憶部６３は、制御部６１の指示に従って、このプログラム６９の実行の際に用いられるデータを制御部６１に供給し、制御部６１から供給されたデータを記憶する。

操作部６４は、キーボード及びマウスなどのポインティングデバイス等と、キーボードおよびポインティングデバイス等を内部バス６０に接続するインターフェイス装置から構成されている。操作部６４を介して、操作者が操作した内容に関する情報が制御部６１に入力される。

表示部６５は、ＣＲＴ（Cathode Ray Tube）またはＬＣＤ（Liquid Crystal Display）などから構成され、操作者が操作情報を入力する場合は、操作用の画面が表示される。

入出力部６６は、データ取得部１０からｋ空間データ１５を入力し、高解像度画像２５を出力するインターフェイスとなる。

データ復元部２０では、制御部６１がプログラム６９を実行することにより、各種機能を実現する。

なお、データ取得部１０において演算を行う部分のハードウエア構成も、図９のデータ復元部２０のハードウエアと同様であってもよいし、共通化されていてもよい。

次に、本発明の実施の形態に係る磁気共鳴画像高速再構成法（磁気共鳴イメージング装置１の動作）について説明する。

図１０に示すように、まず、学習部２２が、教師データを用いた学習を行う（ステップＳ１；深層学習ステップ）。上述のように、学習部２２は、高解像度画像２５のｋ空間データの教師データ５０（図６参照）と、教師データ５０が間引きされて生成されたｋ空間データ１５を用いて、ＭＲＣＮＮ２１を学習させる深層学習を行い、ＭＲＣＮＮ２１の各種パラメータを最適化する。なお、この深層学習ステップでは、所定の空間周波数よりも低周波成分のデータ全体を残し、所定の空間周波数よりも高周波成分についてはランダムにデータを残すアンダーサンプリングフィルタ４０を用いて、教師データ５０が間引きされ、ｋ空間データ１５が生成される。また、この深層学習ステップでは、教師データ５０と、ＭＲＣＮＮ２１を用いて復元された高解像度画像２５に対応するｋ空間データとの誤差に基づいて、誤差逆伝播法による学習が行われる。

続いて、データ取得部１０は、上記シーケンスに従って動作し、磁気共鳴により観測される核磁気共鳴信号（原信号）に基づくｋ空間データ１５を取得する（ステップＳ２；データ取得ステップ）。これにより、被検体の断層画像を得るためのｋ空間データ１５が取得される。なお、このデータ取得においては、図７に示すアンダーサンプリングフィルタ４０に従って、シーケンスの繰り返し数又は観測するスライス位置Ｚの数が規定されている。すなわち、ｋ空間における低周波数成分についてはデータが間引きせずにサンプリングされ、高周波成分についてはランダムに間引きされた状態でデータがサンプリングされる。これにより、データ取得時間（撮像時間）が短縮されている。

続いて、データ復元部２０が、多重解像度畳み込みニューラルネットワーク（ＭＲＣＮＮ２１）を用いて、磁気共鳴により観測される原信号に基づくｋ空間データ１５から高解像度画像２５を復元する（ステップＳ３；データ復元ステップ）。具体的には、ＭＲＣＮＮ２１は、図４に示すように、ｋ空間データ１５から、サイズが異なる（スケールダウンした）ｋ空間データ（１／２データ、１／４データ）を生成する。そして、ＭＲＣＮＮ２１は、ｋ空間データ１５、１／２データ及び１／４データ各々から逆フーリエ変換される画像を、ＳＲＣＮＮ３５を用いて復元する。さらに、ＭＲＣＮＮ２１は、復元されたデータを、アップサンプリング、すなわちデータのサイズを合わせて重ね合わせることにより、高解像度画像２５を復元する。

図１１には、劣化した被検体の断層画像が示され、図１２には、本実施の形態に係る磁気共鳴画像高速再構成法において復元された被検体の断層画像が示されている。図１１及び図１２を比較するとわかるように、被検体の断層画像が高精細に復元されている。なお、この断層画像は、Ｔ１強調画像である。Ｔ１強調画像は、２次元で脳の組織を見るための画像である。Ｔ１強調画像は、大脳皮質又は骨などの解剖学的な構造を見る時に使用される。Ｔ１は、縦磁化の緩和時間である。Ｔ１強調画像は、ＲＦパルスが出射される間隔である繰り返し時間を短くしてＴ１の影響を強調し、ＲＦパルスを出射してからエコー信号の中央までの時間であるエコー時間を短くしてＴ２（横磁化の緩和時間）の影響を少なくすることにより得られる画像である。

本実施の形態に係る磁気共鳴画像高速再構成法による画像の復元状態を評価した。データとしては、頭部のＭＲＩ断層画像（サイズ：２５６×２５６）を用いた。ｔｒａｉｎを１００枚、ｔｅｓｔを５０枚として、本実施の形態に係るＭＲＣＮＮ２１によるデータ復元、ＳＲＣＮＮ３５単独でのデータ復元との間での性能比較を行った。画像の復元の評価方法としては、ＰＳＮＲ、ＭＩＳＳＩＭ、ＫＰＳＮＲが用いられた。

まず、ＰＳＮＲについて説明する。ｘを原画像、ｙを比較したい画像とすると、ＰＳＮＲは、以下の式で定義される。

上式により、ＰＳＮＲの値が大きいほど原画像に近いという評価になる。
ＳＳＩＭは、画像の小領域毎に計算される。計算式を以下に示す。
ただし、
ここで、ｘ，ｙは、それぞれの小領域の画素であり、μ_ｘ，μ_ｙは、それぞれの平均画素値であり、σｘ，σｙは、それぞれの標準偏差であり、σ_ｘｙはｘとｙの共分散である。また、Ｃ_１＝（Ｋ_１Ｌ）^２ _，Ｃ_２＝（Ｋ_２Ｌ）^２であり、Ｌはダイナミックレンジ、Ｋ１、Ｋ２は、任意の値であり、例えば、Ｋ_１＝０．０１、Ｋ_２＝０．０３とすることができる。

ここで、α＝β＝γ＝１かつＣ_３＝Ｃ_２／２とすると、ＳＳＩＭ（ｘ，ｙ）は以下の式で表される。
ＳＳＩＭは小領域ごとに計算され、画像に対するＳＳＩＭの平均を取ったものをＭＳＳＩＭという。ＭＳＳＩＭを次式に表す。

次に、ＫＰＳＮＲについて説明する。ＫＰＳＮＲは周波数領域ごとの復元度を調べるために考案した手法で、Ｋ空間でＰＳＮＲを測定する。図１３、図１４、図１５の白の領域をそれぞれ低周波領域、中間領域、高周波領域として抽出し、それぞれの領域でＫＰＳＮＲを測定した。各領域での測定結果を、ＫＰＳＮＲ（ｌｏｗ）、ＫＰＳＮＲ（ｍｉｄｄｌｅ）、ＫＰＳＮＲ（ｈｉｇｈ）とした。

図１６に示すように、ＭＲＣＮＮ２１は、ＰＳＮＲ、ＭＳＳＩＭ、ＫＰＳＮＲ（ｌｏｗ）、ＫＰＳＮＲ（ｍｉｄｄｌｅ）、ＫＰＳＮＲ（ｈｉｇｈ）全てにおいて、ＳＲＣＮＮ３５単独の場合よりも、評価値が良好であった。すなわち、畳み込みニューラルネットワークの多重化を行った効果、本実施の形態に係るＭＲＣＮＮ２１の効果を確認することができた。

また、本実施の形態に係る磁気共鳴画像高速再構成法を、ＭＲＡ（Magnetic Resonance Angiography）画像に適用した。ＭＲＡ画像とは、３次元の血管の構造を把握するために使用される画像である。ＭＲＡ画像は、ＭＩＰ（Maximum Intensity Projection；最大投影法）処理を行うことにより生成され、血管が強調された画像となる。

ＭＩＰ処理とは、３次元の画像データを任意の方向（投影方向）にある２次元の投影面に投影した場合に、３次元の画像データにおける投影面の同じ位置に投影される複数の画素の最大値を投影面上の各画素の値とする処理のことである。この処理により、３次元の画像データにおいて投影方向に沿った複数のボクセルの中で最も信号強度が高いボクセルの信号強度が、ＭＲＡ画像のピクセルの信号強度となる。

ここでは、患者１人当たりのＭＲＡ画像のデータのサイズを、１７５×５１２×５１２（スライス×縦×横）とした。また、訓練データを２００枚とし、テストデータを１００枚とした。フィルタとしては、図７のアンダーサンプリングフィルタを用いた。また、学習時のハイパーパラメータとして、学習係数（誤差逆伝播の係数）を１．０×１０^−５とし、反復回数を１．０×１０^４とし、バッチサイズを１０とした。

図１７（Ａ）〜図１７（Ｃ）には、頭頂から見た場合の、すなわち頭頂から見た場合に視線方向に投影されるＭＩＰ画像の一例が示されている。図１７（Ａ）は、元画像（正解画像）である。この正解画像は、時間をかけて計測された高精細な画像である。ここで、この正解画像を教師データ５０としてＭＲＣＮＮ２１により学習を行った。例えば、アンダーサンプリングフィルタ４０を用いて図１７（Ａ）に示す元画像から、図１７（Ｃ）に示す劣化画像を生成し、ＭＲＣＮＮ２１の学習を行った。さらに、学習されたＭＲＣＮＮ２１により、図１７（Ｃ）に示す劣化画像から、図１７（Ｂ）に示す復元画像が得られた。

図１８（Ａ）〜（Ｃ）には、図１７（Ａ）〜（Ｃ）の白枠部分の拡大図が示されている。図１８（Ａ）と図１８（Ｂ）とを比較して示すように、ＭＲＣＮＮ２１で復元された復元画像は、元画像（正解画像）と比べ、精細さに遜色はなかった。

図１９（Ａ）〜図１９（Ｃ）には、横方向から見た場合のＭＩＰ画像の一例が示されている。図１９（Ａ）は、元画像（正解画像）である。この正解画像は、時間をかけて高精細な画像である。同様に、この正解画像を教師データ５０としてＭＲＣＮＮ２１により学習を行った。例えば、アンダーサンプリングフィルタ４０を用いて図１７（Ａ）に示す元画像から、図１９（Ｃ）に示す劣化画像を生成し、ＭＲＣＮＮ２１の学習を行った。そして、学習されたＭＲＣＮＮ２１により、図１９（Ｃ）に示す劣化画像から、図１９（Ｂ）に示す復元画像が得られた。

図２０（Ａ）〜（Ｃ）には、図１９（Ａ）〜（Ｃ）の白枠部分の拡大図が示されている。図２０（Ａ）と図２０（Ｂ）とを比較して示すように、ＭＲＣＮＮ２１で復元された復元画像は、元画像（正解画像）と比べ、精細さに遜色はなかった。

図２１（Ａ）〜図２１（Ｃ）には、正面方向から見た場合のＭＩＰ画像の一例が示されている。図２１（Ａ）は、元画像（正解画像）である。この正解画像は、時間をかけて高精細な画像である。ここで、この正解画像を教師データ５０としてＭＲＣＮＮ２１により学習を行った。例えば、アンダーサンプリングフィルタ４０を用いて図２１（Ａ）に示す元画像から、図２１（Ｃ）に示す劣化画像を生成し、ＭＲＣＮＮ２１の学習を行った。そして、学習されたＭＲＣＮＮ２１により、図２１（Ｃ）に示す劣化画像から、図２１（Ｂ）に示す復元画像が得られた。

図２２（Ａ）〜（Ｃ）には、図２１（Ａ）〜（Ｃ）の白枠部分の拡大図が示されている。図２２（Ａ）と図２２（Ｂ）とを比較して示すように、ＭＲＣＮＮ２１で復元された復元画像は、元画像（正解画像）と比べ、精細さに遜色はなかった。

なお、画像の復元度の評価を行うために、上記式（１４）のＰＳＮＲが使用された。図２３には、ＰＳＮＲの結果が示されている。具体的には、ＭＩＰ処理を行う前のＭＲＡ原原画像について、元画像に対する劣化画像のＰＳＮＲ値と、元画像に対する復元画像のＰＳＮＲ値と、復元画像のＰＳＮＲ値とが求められ、それらの差が復元度として求められた。ＭＩＰ画像（頭頂）と、ＭＩＰ画像（横方向）と、ＭＩＰ画像（正面方向）とについても、同様にＰＳＮＲ値が求められた。図２３に示すように、ＭＩＰ画像（頭頂）と、ＭＩＰ画像（横方向）と、ＭＩＰ画像（正面方向）とについて、それぞれの復元度は大幅に向上しており、ＭＩＰ処理前のＭＲＡ原画像と、ＭＩＰ処理後の画像とを比較すると、ＭＩＰ処理を行った方がより復元度が向上していることが分かった。これは、ＭＩＰ処理後の画像は、投影方向に沿ったボクセルの最大値に基づく画像であり、ＭＩＰ処理前に比べノイズ成分が低減されているためであると考えられる。

以上詳細に説明したように、本実施の形態によれば、処理時間が圧縮センシングよりも短い多重解像度畳み込みニューラルネットワーク（ＭＲＣＮＮ２１）を用いて観測点数の少ない第１のｋ空間データ１５をスケールダウンし、多重化して超解像処理を行って得られる画像を用いて画像の再構成を行う。このため、観測点数を少なくして撮像時間を短縮しつつ、多重解像度畳込みニューラルネットワークを用いて再構成時間を短縮して、高周波成分だけでなく低周波成分も損なうことなく高解像度画像２５を復元することができる。この結果、撮像時間及び再構成時間を双方短縮しつつ、高解像度画像２５を復元することができる。

ＭＲＩ断層撮影では、周波数領域でアンダーサンプリングが行われるため、様々な周波数成分が不足してしまう。また、ＳＲＣＮＮ３５単体では、フィルタのサイズが小さいため、高周波成分しかとらえることができないという問題があり、フィルタを大きくするしかなく、この場合にはパラメータ数が増えて学習が困難になる。そこで、本実施の形態に係る磁気共鳴イメージング装置１では、超解像畳み込みニューラルネットワークの多重化を行って、入力画像のサイズを小さくすることで相対的にフィルタのサイズを大きくしたので、低周波成分の復元も行うことができるようになり、この結果、高精細な画像を復元することができるようになった。

また、本実施の形態では、画像ではなく、磁気共鳴信号（原信号）から構成されるｋ空間そのものに対する学習により、データ復元を行う。これにより、低域から高域までの様々な空間周波数成分を直接復元することができるようになるので、より高精細な画像の復元が可能となる。このため、これまで発見が困難であった病変の発見も可能になることが期待される。

また、本実施の形態では、画像の復元に有益な低周波成分全体を残し、病変を発見する上で重要な高周波成分をランダムに残すアンダーサンプリングフィルタを用いて学習を行った。これにより、病変を発見し易く、かつ、高精細な画像を復元することが可能となった。

なお、本実施の形態では、アンダーサンプリングフィルタ４０のデータ取得率を１０％としたが、本発明はこれには限られない。データ取得率は、１０％を超えてもよいし、１０％未満であってもよい。

なお、アンダーサンプリングフィルタ４０の形状は、図７に示すものには限られない。アンダーサンプリングフィルタ４０の形状は、原信号の低周波成分と、高周波成分とが大きく損なわれないようなものであるのが望ましい。

また、学習部２２は、高解像度ｋ空間データ１５と教師データ５０との誤差に基づく誤差逆伝播法により学習を行った。これにより、ｋ空間データにおける誤差に基づいて直接学習が行われるようになるので、学習時間を短縮することが可能となる。なお、学習方法は、誤差逆伝播法には限られない。ＭＲＣＮＮ２１のパラメータを最適化できるのであれば、様々な学習方法を適用可能である。

なお、上記実施の形態では、ＭＲＣＮＮ２１における多重化を３つとしたが、２つでもよいし、４つ以上でもよい。多重化が２つの場合には、１／２データが生成され、多重化が４つ以上の場合には、１／２データ、１／４データ、１／８データ、・・・を生成して、それぞれをＳＲＣＮＮ３５に入力すればよい。

また、多重化するのはＳＲＣＮＮ３５に限られない。他の畳み込みニューラルネットワークを多重化するようにしてもよい。

なお、本発明は、データ取得（撮影）のパルスシーケンスに制限はない。スピンエコー法、高速スピンエコー法、グラディエントエコー法、エコープレナー法、ハーンエコー法、スティミュレイテドエコー法など、様々な方法に対応することができる。

その他、データ復元部２０のハードウエア構成やソフトウエア構成は一例であり、任意に変更および修正が可能である。

制御部６１、主記憶部６２、外部記憶部６３、操作部６４、表示部６５及び入出力部６６、内部バス６０などから構成されるデータ復元部２０の処理を行う中心となる部分は、上述のように、専用のシステムによらず、通常のコンピュータシステムを用いて実現可能である。例えば、前記の動作を実行するためのコンピュータプログラムを、コンピュータが読み取り可能な記録媒体（フレキシブルディスク、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ等）に格納して配布し、当該コンピュータプログラムをコンピュータにインストールすることにより、前記の処理を実行するデータ復元部２０を構成してもよい。また、インターネット等の通信ネットワーク上のサーバ装置が有する記憶装置に当該コンピュータプログラムを格納しておき、通常のコンピュータシステムがダウンロード等することで磁気共鳴イメージング装置１を構成してもよい。

コンピュータの機能を、ＯＳ（オペレーティングシステム）とアプリケーションプログラムの分担、またはＯＳとアプリケーションプログラムとの協働により実現する場合などには、アプリケーションプログラム部分のみを記録媒体や記憶装置に格納してもよい。

搬送波にコンピュータプログラムを重畳し、通信ネットワークを介して配信することも可能である。たとえば、通信ネットワーク上の掲示板（BBS, Bulletin Board System）にコンピュータプログラムを掲示し、ネットワークを介してコンピュータプログラムを配信してもよい。そして、このコンピュータプログラムを起動し、ＯＳの制御下で、他のアプリケーションプログラムと同様に実行することにより、前記の処理を実行できるように構成してもよい。

この発明は、この発明の広義の精神と範囲を逸脱することなく、様々な実施の形態及び変形が可能とされるものである。また、上述した実施の形態は、この発明を説明するためのものであり、この発明の範囲を限定するものではない。すなわち、この発明の範囲は、実施の形態ではなく、特許請求の範囲によって示される。そして、特許請求の範囲内及びそれと同等の発明の意義の範囲内で施される様々な変形が、この発明の範囲内とみなされる。

本発明は、磁気共鳴画像の撮影時間、データ処理時間の短縮及び高解像度化に適用することができる。

１磁気共鳴イメージング装置、１０データ取得部、１５ｋ空間データ、２０データ復元部、２１多重解像度畳み込みニューラルネットワーク（ＭＲＣＮＮ）、２２学習部、２５高解像度画像、３５超解像畳み込みニューラルネットワーク（ＳＲＣＮＮ）、４０アンダーサンプリングフィルタ、４１フーリエ変換部、４５パラメータ設定部、５０教師データ、６０内部バス、６１制御部、６２主記憶部、６３外部記憶部、６４操作部、６５表示部、６６入出力部、６９プログラム

Claims

多重解像度畳み込みニューラルネットワークを用いて、磁気共鳴により観測される原信号に基づく第１のｋ空間データから高解像度画像を復元するデータ復元ステップを含み、
前記多重解像度畳み込みニューラルネットワークでは、
前記第１のｋ空間データから、サイズが異なる第２のｋ空間データを生成し、
前記第１のｋ空間データ及び前記第２のｋ空間データ各々から得られる画像を、超解像畳み込みニューラルネットワークを用いて超解像化し、
超解像化された画像を、サイズを合わせて重ね合わせることにより、前記高解像度画像を復元する、
磁気共鳴画像高速再構成法。
前記データ復元ステップに先立って、
高精細なｋ空間データを教師データとし、前記教師データが間引きされて生成されたｋ空間データを前記第１のｋ空間データとして、前記多重解像度畳み込みニューラルネットワークを学習させる深層学習ステップを含む、
請求項１に記載の磁気共鳴画像高速再構成法。
前記深層学習ステップでは、
所定の空間周波数よりも低周波成分のデータ全体を残し、前記空間周波数よりも高周波成分についてはランダムにデータを残すアンダーサンプリングフィルタを用いて、前記教師データを間引きすることにより、学習用の前記第１のｋ空間データを生成する、
請求項２に記載の磁気共鳴画像高速再構成法。
前記深層学習ステップでは、
前記教師データと、前記多重解像度畳み込みニューラルネットワークを用いて復元された前記高解像度画像に対応する第３のｋ空間データとの誤差に基づく誤差逆伝播法を用いて学習を行う、
請求項２又は３に記載の磁気共鳴画像高速再構成法。
磁気共鳴により観測される原信号に基づく第１のｋ空間データを取得するデータ取得部と、
多重解像度畳み込みニューラルネットワークを用いて、前記第１のｋ空間データから高解像度画像を復元するデータ復元部と、
を備え、
前記多重解像度畳み込みニューラルネットワークは、
前記第１のｋ空間データから、サイズが異なる第２のｋ空間データを生成し、
前記第１のｋ空間データ及び前記第２のｋ空間データ各々から得られる画像を、超解像畳み込みニューラルネットワークを用いて超解像化し、
超解像化された画像を、サイズを合わせて重ね合わせることにより、前記高解像度画像を復元するように構成されている、
磁気共鳴イメージング装置。