エッジ情報と LiGME モデルを用いたマルチコントラスト圧縮センシング MRIMulti-Contrast Compressed Sensing MRI Using Edge Information with LiGME Model
北原 大地

加藤 里佳子黒田 大貴
立命館大学 情報理工学部

平林 晃

Daichi KITAHARA
Rikako KATO
Hiroki KURODA
Akira HIRABAYASHI
College of Information Science and Engineering, Ritsumeikan University
あらまし

体内を非侵襲で撮像する方法の 1 つとして MRI

不鮮明な画像が再構成されてしまう．削減されたデータがある．MRI は検査に時間がかかるため，取得データ数をからも鮮明な画像を再構成するために，圧縮センシング

減らして撮像時間を短縮した上で，なるべく高精度な画像(Compressed Sensing: CS) 理論が利用されている [4]–[8]．

再構成を目指す，圧縮センシング MRI が研究されている．

圧縮センシング [9] とは，未知の対象信号が適切な表現

臨床現場では一度の検査で複数種類の MR 画像を撮像する空間においてスパースになることを利用して，少数の観測

が，Ehrhardt らは 1 種類の MR 画像のみ全データを取得し信号から対象信号を高精度に復元する技術である．このてそのエッジ情報を活用する，マルチコントラスト圧縮技術を利用し，Lustig らは圧縮センシング MRI (CSMRI)

センシング MRI を提案した．他の種類の画像は，共通のと呼ばれる高速 MR 撮像法を提案した [4]．この手法では，

エッジ情報に基づき少数の観測データから高精度に再構成

観測信号との整合性，離散ウェーブレット係数のスパースされる．本研究では，Ehrhardt らの手法を近年提案された性，及び差分値のスパース性を評価する損失関数を定義し，

LiGME モデルに変更することを提案する．LiGME モデル

この値を最小化することで MR 画像を再構成する．圧縮率を用いることで，エッジ情報を考慮した全変動がグループ

1/3 程度の観測信号からでも鮮明な画像が再構成される．

ℓ0 擬ノルムに近づくため，大きなエッジを復元しやすく

他にも，学習用画像や推定途中の対象画像から，MR 画像のなる．実データを用いた実験で，提案法の有効性を示す． スパース表現自体も求める手法が提案されている [5]–[7]．
Ehrhardt らは，検査時に同一患者の同一部位を複数の
1 はじめにコントラストで画像化することに着目して，1 種類の画像
磁気共鳴画像法 (Magnetic Resonance Imaging: MRI) [1], のみにフルサンプリングを行い，他のコントラストの画像[2] は強力な磁力を利用して生体の内部構造を画像化するには間引いたサンプリングを行って検査時間全体を短縮装置である．コンピュータ断層撮影法 (Computed Tomog-

する，マルチコントラスト CSMRI を提案した [8]．この

raphy: CT) [3] と比較して，X 線を被曝しないという利点

手法では，フルサンプリングの観測値から得られた高品質がある一方で，検査にかかる時間が長いという欠点がある． な MR 画像のエッジ情報を用いて，指向性全変動と呼ば臨床現場における MRI 検査では，放射線科医が適切な れる正則化項を定義する．そして，観測信号との二乗誤差
診断を下すのに十分な情報を得るために，一度の検査でと指向性全変動の和を最小化することで，異なるコント

繰り返し時間 (TR) とエコー時間 (TE) の設定を変えて複数

ラストの画像が少数の観測値から高精度に再構成される．

種類の MR 画像を取得している．取得頻度が高い画像と近年，ℓ1 ノルムに代わるスパース性の評価尺度として，して，解剖学的構造の描出に長けた T1 強調画像と，病変Generalized Minimax Concave (GMC) 関数が Selesnick にの描出に長けた T2 強調画像があり (第 2.3 節の図 1 参照)， より提案された [10]．GMC 関数は ℓ0 擬ノルムと似た形状
どちらも撮像には 2 分から 3 分程度の時間を必要とする． を持つ非凸な連続関数であり，ℓ1 ノルムよりもスパース性
CT は数十秒で撮像を終えることから，患者の負担を考慮を正確に評価できる．安倍らは，GMC 関数を一般化した
して，MRI の撮像時間を短縮することが求められている． Generalized Moreau Enhanced (GME) 関数を定義し，観測MRI では，観測データを時間軸上で順次取得するため， 信号との二乗誤差と線形変換を組み込んだ GME 関数の和取得データ数を減らすことで撮像時間を短縮できる．削減の厳密な最小化を行う，Linearly Involved GME (LiGME)

された観測データに対して通常の画像再構成手法である

モデルを提案し，その優れた信号復元性能を示した [11]．

逆フーリエ変換を適用すると，アーチファクトを含んだ

本研究では，Ehrhardt らの手法 [8] が LiGME モデルに本研究の一部は JSPS 科研費 (19K20361) の御支援を受けて行われた．

変更可能なことを明らかにして，文献 [11] のアルゴリズム

に基づいた MR 画像再構成手法を提案する．LiGME モデル

2 次元離散フーリエ変換行列であり，εfull ∈ CN はモデルを用いることで，指向性全変動がグループ ℓ0 擬ノルムに誤差と観測誤差の和である．雑音のエネルギー ∥εfull ∥22 が近い値で計算されるようになるため，大きなエッジを保存

小さい場合は，観測信号 yfull を逆離散フーリエ変換した

した状態で画像を再構成できると考えられる．実データ

結果を非負実数値 N 次元ベクトルの集合 RN+ に射影するを用いた計算機シミュレーションにより，提案法が従来法ことで，即ち PRN(F H yfull ) によって，MR 画像 x が高精度+と比べて高速かつ高精度に画像を再構成することを示す． に再構成される3 ．ここで，RN: CN → RN+ への射影 PRN++

は，複素数値ベクトルの虚部を切り捨てて，残った実部

従来の MR 画像再構成

2

実数全体の集合を R で，非負実数全体の集合を R+ で，の負値を 0 に置き換えるという処理で簡単に実行できる．

複素数全体の集合を C で表し，虚数単位を ı ∈ C で表す， 2.2

圧縮センシング MRI (CSMRI)

即ち ı = −1 である．複素数 c ∈ C の複素共役を c̄ で表す．

MRI では，観測データ y[k1 , k2 ] が一括ではなく時間軸ベクトルと行列をそれぞれボールド体の小文字と通常書体

上で順次取得される．このため，取得データ数を減らせば2の大文字で表す．N 次の単位行列を IN で表す．ベクトル

撮像時間を短縮できる．M (< N ) 個の観測データ y[k1 , k2 ]及び行列の転置とエルミート転置を (·)T と (·)H で表す． からなる信号を y ∈ CM で表せば，CSMRI での観測式はN 次のエルミート行列 A ∈ CN ×N の正定値性と半正定値
y = Fu x + ε(3)性を A ≻ O と A ⪰ O で表す．N 次元ベクトル x :=NN

(x1 , x2 , . . . , xN ) ∈ √C の ℓ2 ノルム ∥·∥2 と ℓ1 ノルム ∥·∥1∑N∑N
2i=1 |xi | と ∥x∥1 :=i=1 |xi | のTNをそれぞれ ∥x∥2 :=ように定義する．また，GQ = N を満たす 2 つの正整数のようにモデル化される．ここで，Fu ∈ CM ×N は間引き
フーリエ変換行列であり，k 空間のサンプリング点を表す∑N
行列 S := (si,j ) ∈ {0, 1}M ×N (s.t. ∀i j=1 si,j = 1 and∑M∀j i=1 si,j ≤ 1) を用いて，Fu = SF と表される．また，G,Q
G と Q を用いて，グループℓ1 ノルム ∥·∥G,Q
1,2 を ∥x∥1,2 :=∑G∑G √∑Q−12
i=1 ∥xgi ∥2 :=i=1
j=0 |xi+jG | のように定義する． ε ∈ CM は雑音である．式 (3) における観測信号 y の圧縮率

2.1

′
′
をM
N と定義する．取得していないフーリエ係数 y[k1 , k2 ] の

通常の MRIMRI の撮像方式にはスピンエコー (Spin Echo: SE) 系とグラディエントエコー (Gradient Echo: GRE) 系の 2 種類がある．SE 系の場合は MR 画像が非負実数値信号となり，

GRE 系の場合は MR 画像が複素数値信号となる．本研究では SE 系の撮像方式を想定し，N 個の非負実数値で構成される N 次元ベクトル x ∈ RN
+ で MR 画像を表現する．

k 空間 における最も基本的なサンプリングでは，
11y[k1 , k2 ] = √
N

N1
2 −1

N2
2 −1

N

N

∑

∑

x[n1 , n2 ]e−ı2π

部分に 0 や y[−k1′ , −k2′ ] を埋め込んだ信号を ypad ∈ CN で

(F H ypad ) が最も単純な再構成画像となるが4 ，
表せば，PRN+
画像全体が不鮮明であり，折り返しパターンなどのアーチ
ファクトも発生するため，診断に使用することはできない．
この問題を解決するために，Lustig らは凸最適化問題

minimize
x∈RN
+

(k n

k2 n2
1 1
N1 + N2)n1 =− 21 n2 =− 22(1)のように格子点上の複素フーリエ係数 y[k1 , k2 ] ∈ C を観測する2 ．ここで，N1 と N2 は偶数で，それぞれ MR 画像 x の1∥Fu x − y∥22 + λ∥Φx∥1 + µTV(x)2(4)

を解くことで MR 画像 x を再構成することを提案した [4]．
式 (4) の第 1 項はデータ整合性で，第 2 項と第 3 項は MR
画像らしさを評価する正則化項である．SE 系の場合は xが非負実数値信号となるため，通常の自然画像処理と同じ
正則化項を利用できる．Φ ∈ RN ×N は離散ウェーブレット

縦と横の画素数を表し，N1 N2 = N を満たす．また，k1 =変換，TV : RN → R+ は等方性全変動 (Isotropic Total Vari-− N21 , − N21 +1, . . . , N21 −1，k2 = − N22 , − N22 +1, . . . , N22 −1
である．式 (1) の観測モデルは積分の近似式 (数値積分) で

ation) であり，2 次元差分行列
D := [DhT , DvT ]T を用いて5∑N √ 2TV(x) := n=1 ∆xh,n + ∆x2v,n = ∥Dx∥N,2
1,2 と表される．

あり，実際の観測データ y[k1 , k2 ] にはモデル誤差と観測

λ > 0 と µ > 0 は評価関数内の各項の重みを決定するパラ誤差が含まれている．全ての観測データ y[k1 , k2 ] をまとメータである．式 (4) に基づく手法は，圧縮率が 1/3 程度

めたベクトルを yfull ∈ C

までならば SE 系の MR 画像 x を高精度に再構成できる．

Nで表せば，式 (1) の観測式はyfull = F x + εfull(2)のように表現される．ここで，F ∈ CN ×N は正規化された
1 MRI では周波数空間を「k 空間」と呼び，ここでサンプリングを行う．
2 SE 系の場合は，y[0, 0], y[0, − N2 ], y[− N1 , 0], y[− N1 , − N2 ] ∈ R
2
2
2
2となり，k1 ̸= − N21 と k2 ̸= − N22 に対して y[k1 , k2 ] = y[−k1 , −k2 ]，y[− N21 , k2 ] = y[− N21 , −k2 ]，y[k1 , − N22 ] = y[−k1 , − N22 ] も成り立つ．

3 GRE 系の場合には射影計算が不要となり，逆離散フーリエ変換した
結果 F H yfull が再構成画像となる．また，SE 系の場合は，複素フーリエ係数 y[k1 , k2 ] の複素共役の関係 (脚注 2 参照) を利用すれば，式 (2) の約半分の観測信号 yhalf からでも MR 画像 x を高精度に再構成できる．
4 GRE 系の場合は脚注 2 の関係が使えないため，0 を埋め込んだ後に逆離散フーリエ変換した結果 F H ypad = FuH y が単純再構成画像となる．
5 D ∈ {−1, 0, 1}N ×N と D ∈ {−1, 0, 1}N ×N は，それぞれ水平v
h
方向と垂直方向の差分行列を表しており，右端の画素値と左端の画素値の差分，及び，下端の画素値と上端の画素値の差分も巡回して計算する．

のように定義される．ここで，式 (6) の η > 0 は，数値的不安定性を回避するための小さな正定数である．式 (7) の√ξen := 2 − ∥ξn ∥2 ξn は，第 n 画素における参照画像 u の2

重み付きエッジ情報であり，∆un := (∆uh,n , ∆uv,n )T と平行で 0 ≤ ∥ξen ∥2 < 1 となる．∥∆xn ∥2 = ∆x2 + ∆x2
2

2

h,n

v,n

を一定にしたとき，∆xn と ξen が平行ならば dTV(x) の値は最小となり，∆xn と ξen が直交するならば dTV(x)の値は最大となる．式 (5) に基づく手法は，少ない観測値(b) T2 強調画像からでも MR 画像の細部構造を高精度に再構成できる [8]．

図 1: T1 強調画像と T2 強調画像の類似性

3 エッジ情報と LiGME モデルを用いた MR 画像再構成(a) T1 強調画像

2.3

3.1 LiGME モデルに基づく信号復元

エッジ情報を用いたマルチコントラスト CSMRI

実際の臨床において MR 画像を 1 枚のみ撮像することはまずなく，必ず同一部位を異なる設定で複数回撮像する．
これにより，同一スライスでコントラストが異なる複数のMR 画像が得られ，医師はそれらを見比べることで診断を行う．Ehrhardt らはこの診察過程に着目し，ある設定では式 (2) のフルサンプリングを，他の設定では式 (3) の間引き
サンプリングを行うことで，画質を保ちつつ診察時間全体を短縮するマルチコントラスト CSMRI を提案した [8]．
ここでは，簡単のため，T1 強調画像をフルサンプリング
で撮像し，T2 強調画像を間引きサンプリングで撮像するものとする．雑音のエネルギーが小さい場合には，T1 強調
ref
画像がフルサンプリングの観測値から u := PRN(F H yfull)+のように簡単に再構成される．T1 強調画像と T2 強調画像
では，コントラストが大きく異なるもののエッジの構造は共通であるため (図 1)，撮像された高品質な T1 強調画像u ∈ RN+ を参照画像として利用することで，T2 強調画像
Mから高精度に再構成
x ∈ RN+ が間引いた観測値 y ∈ C

可能となる．具体的に，Ehrhardt らは以下の凸最適化問題

1∥Fu x − y∥22 + µ dTV(x)2minimize
x∈RN
+

(5)

を解くことで，MR 画像 x を再構成した．ここで，dTV :R

N→ R+ は指向性全変動 (Directional Total Variation) で，

参照画像 u の第 n 画素での水平・垂直方向の正規化差分値
[
]
[
]
∆uh,n
ξn,1
1
(6)
ξn :=
:= √
ξn,2
∆u2h,n + ∆u2v,n + η 2 ∆uv,nを用いて，

dTV(x) :=

N
∑
n=1

=

N
∑

(I2 − ξn ξnT )
[

n=1

=

N √
∑

2
1 − ξn,1−ξn,1 ξn,2

−ξn,1 ξn,22
1 − ξn,2

2

][

∆xh,n
∆xv,n

∆x2h,n + ∆x2v,n − (ξenT ∆xn )2n=1を解くことで観測信号 y ∈ RM から未知信号 x ∈ RN を推定する．ここで，L ∈ RK×N はある行列，Ψ : RK → R
は下半連続なある真凸関数である．圧縮センシング理論[9] に基づく場合は，Ψ が ℓ1 ノルムやその類似関数となる．
圧縮センシング理論の目的は，少数の観測データからスパースな未知信号を推定することであるが，スパース性を正確に評価するためには，ℓ1 ノルムよりも ℓ0 擬ノルム
や ℓp 準ノルム (0 < p < 1) を用いた方がよい．しかし，Ψ
に ℓ0 擬ノルムや ℓp 準ノルムを用いると，式 (8) の問題は非凸になり，大域的最適解を求めることが困難となる．

Selesnick は，L が単位行列 IN で Ψ が ℓ1 ノルムのとき，][12ΨB (w) := Ψ(w) − min Ψ(v) + ∥B(w − v)∥2 (9)2v∈RKによって定義される Generalized Minimax Concave (GMC)
関数を，ℓ1 ノルムに代わる正則化項として式 (8) の問題で利用することを提案した [10]．GMC 関数は ℓ0 擬ノルムに近い形状を持つ連続関数であり，B ∈ RJ×K はその形状を決定する行列である．GMC 関数自体は非凸であるものの，

AT A − µB T B が半正定値行列であるならば，最適化問題
全体は凸になり，大域的最適解を求めることが可能となる．
安倍らは，式 (9) を強圧的かつ対称性 Ψ(−w) = Ψ(w)を満たす全ての凸関数 Ψ に拡張して，これを Generalized

Moreau Enhanced (GME) 関数と呼んだ．更に，任意の行列
L も組み込まれた Linearly Involved GME (LiGME) モデル6
1
minimize ∥Ax − y∥22 + µΨB (Lx)
(10)
N
2
x∈R

[
]
∆xh,n
∆xv,n

凸最適化に基づく信号復元手法の多くが，以下の問題
1
minimize ∥Ax − y∥22 + µΨ(Lx)
(8)
N
2
x∈R

を提案した．式 (10) の問題が凸になるための十分条件は]AT A − µLT B T BL ⪰ ON

(11)

2

であり，大域的最適解を求めることにも成功している [11]．

(7)

1
2
B (Lx) = Ψ(Lx)−minv [Ψ(v)+ 2 ∥B(Lx−v)∥2 ]
1と似た (Ψ ◦ L)B (x) = Ψ(Lx) − minv [Ψ(Lv) + 2 ∥B(x − v)∥22 ] も提案されているが，A = IN で Ψ◦L も 1 次元全変動のみに限定される．
6 文献 [12] で，Ψ

3.2

LiGME モデルによるマルチコントラスト CSMRIや√本研究では，式 (5) の Ehrhardt らの最適化問題を LiGME
モデルに変更することで，再構成精度の更なる向上を図る．
まず，式 (5) の問題が式 (8) の形式で表現できることを示す．
式 (7) の dTV は，全ての n で ξen ̸= 0 が成り立つとき，dTV(x)
=

N
∑

√

e 2
 1 − ∥ξn ∥2

n=1



ξen,1
∥
  ξen ∥2e
1 − ξen,2∥ ξn ∥ 2



ξen,2
∆xh,n
∥ξen ∥2  

ξen,1
∆x
v,n
∥ξen ∥2

B=[b†F
θ SH1 Dh
2µ]
RT W −1 ∈ C2M ×2N

b v† FSH2 D(17)などが挙げられる．式 (16) と式 (17) において，θ ∈ (0, 1]
b † ∈ CN ×N と D
b † ∈ CN ×N はそれぞれであり，対角行列 Dhv

bh と Db v の疑似逆行列である．また，H1 ∈ [0, 1]N ×N とD

H2 ∈ [0, 1]N ×N は周波数応答を意味する対角行列であり，2

文献 [12] ではハイパスフィルタにすることを勧めている．

(12)と表される．式 (12) から dTV は，水平・垂直方向の差分値を計算した後に，各画素で差分ベクトルを回転させてξen に平行な成分と直交する成分に分解し，平行な成分に√小さな重み 1 − ∥ξen ∥22 ∈ (0, 1] を乗じてから ℓ2 ノルムを計算していることが分かる．なお，ξen = 0 となる n が存在しても，そのような n に対しては差分ベクトル ∆xn の ℓ2

3.3

画像再構成アルゴリズム

提案する式 (13) の凸最適化問題には，式 (10) の LiGME
モデルと異なる点が存在する．式 (13) の問題では，x ∈ RN+という制約条件が新たに付与されているのである．実は，
式 (13) のような制約付き LiGME モデルに対しても，制約集合への射影が計算可能ならば，文献 [11] のアルゴリズム

ノルムをそのまま計算するだけなので，問題にならない． を僅かに修正するだけで，最適解を求めることができる．
g : RN → (−∞, ∞] を近接写像が計算可能な凸関数とし，以上から，水平方向と垂直方向の差分行列 Dh ∈ RN ×N と1
2
Dv ∈ RN ×N ，回転行列に相当する直交行列 R ∈ R2N ×2N ， 2 ∥Ax − y∥2 + µΨB (Lx) + g(x) の最小化問題を考える．
重み行列に相当する対角行列 W ∈ (0, 1]2N ×2N を用いて， このとき，∀x Ψ(Lx) < ∞ または ∀x g(x) < ∞ であれば，
[ ]
(
[ ] )
文献 [11] のアルゴリズムで x の更新式の最後に prox 1 g の
N,2
σ
Dh
Dh
計算を追加するだけで，最適解を厳密に求めることが可能
dTV(x) = W R
x=: Ψ W R
x
Dv
Dv
1,2
となる7 ．式 (13) の凸最適化問題は，A := Fu ∈ CM ×N ，Nのように dTV を表すことができる．式 (5) の最適化問題は L := W RD ∈ R2N ×N ，Ψ := ∥·∥N,2
1,2 で，g を R+ の指示([ ] )関数とした場合に一致するので，再構成アルゴリズムはDh1minimize ∥Fu x − y∥22 + µΨ W R
x

N
2
x∈R+
Dv
x
[xk − σ1 (AH A − µLT B H BL)xk
 k+1 = PRN+となるため，これを LiGME モデルに変更した問題は

− σµ LT B H Bvk − σµ LT wk + σ1 AH y]
(
[ ] )

Dh
1

x
(13)  vk+1 = prox µ Ψ [ µτ B H BL(2xk+1 − xk ) + vk − µτ B H Bvk ]minimize ∥Fu x − y∥22 + µΨB W R
N
τ

2
x∈R+
Dv
wk+1 = proxΨ∗ [L(2xk+1 − xk ) + wk ]となる．式 (13) が本研究で提案する最適化問題である．
(18)
提案する式 (13) の問題の凸性を保証するには，式 (11) でのようになる．ここで，パラメータ σ > 0 と τ > 0 は，あるA := Fu = SF かつ L := W RD = W R[DhT , DvT ]T としたκ > 1 に対して，それぞれ σIN − κ2 AH A − µLT L ≻ ON
[ ]
]
[
Dh
と τ I2N − µ( κ2 + κ2 )B H B ⪰ O2N を満足する必要がある．
F H S T SF − µ DhT DvT RT W B H BW R
⪰ ON
Dv(14) 4 計算機シミュレーションが成立するように式 (9) の行列 B を設計する必要がある．

文献 [8] の実験で使われた実際の頭部画像「patient B」とDh と Dv は Block Circulant with Circulant Blocks (BCCB)

サンプリングマスク「cartesianX random 0 32」を本研究

行列であり，離散フーリエ変換 F を用いて対角化できるでも用いた．頭部の T1 強調画像と T2 強調画像は第 2.3 節

b h ∈ CN ×N と D
b v ∈ CN ×N を用いてため，複素対角行列 Dの図 1(a) と図 1(b) に示されており，どちらも画像サイズは

b hF
Dh = F H D

かつbvF
Dv = F H D(15)のように表現することができる [5]．式 (15) の表現を利用すれば，式 (14) の条件を満たす行列 B として，√]θ [b † F SH2 D
b v† F RT W −1 ∈ CM ×2N
B=
SH1 D
h
4µ
(16)

242×242 = 58,564 = N である．ここで，画素値を 255 で割ることで信号値を 1 以下にした．サンプリングマスク
は図 2(a) に示されており，取得データ数が M = 8, 480 であるため，圧縮率は M
N = 0.1448，即ち 1/7 程度となる．
提案法の有効性を示すため，従来法である Ehrhardt らの7 文献 [11] と文献 [13], [14] の凸最適化問題とアルゴリズム及び最適解への収束の証明を組み合わせることで，比較的簡単に確かめられる．

表 1: T1 強調画像を参照し T2 強調画像を再構成した結果
検証項目 / 手法 初期画像
従来法
提案法

PSNR

22.81 dB

26.78 dB

28.09 dB

再構成時間—

7.85 秒

4.22 秒

表 2: T2 強調画像を参照し T1 強調画像を再構成した結果
検証項目 / 手法 初期画像
従来法
提案法(a) サンプリングマスク(b) 部分フーリエ法適用後

図 2: SE 系を想定した場合の k 空間のサンプリングマスク

PSNR

21.36 dB

23.63 dB

24.43 dB

再構成時間—

8.48 秒

4.12 秒マルチコントラスト CSMRI を LiGME モデルに変更する

手法 [8] と再構成精度及び再構成時間を比較する．実験

ことで，エッジの向きと大きさをより正確に復元できる

は，m-Book T510XN-M2SH5 (Windows 10, Intel Core i7-ようになる．実際の頭部画像を用いた数値実験で，従来法8750H, 2.20 GHz, 32 GB) 上で MATLAB R2019a を用いてと比べて，提案法では精度と計算時間がともに改善した．

行われた．従来法には文献 [8] のプログラムを使用した．

参考文献

まず，T1 強調画像を参照画像 u として，T2 強調画像を間引かれた観測信号 y から再構成した．図 2(b) のように，フーリエ係数の複素共役の関係に基づいて値を埋め込んだ(部分フーリエ法と呼ばれる) 信号 ypad から，逆フーリエ変換によって初期画像 PRN(F H ypad ) を作成し，従来法と+提案法のアルゴリズムで画像の更新を行った．文献 [8] のプログラムによって，従来法の各パラメータは良好に設定されている．提案法では，µ = 0.003，η = 0.035 とした上で，行列 B には式 (17) を用いて，θ = min(diag(W 2 ))，H1 = H2 = IN と設定した．式 (18) のパラメータ σ と τは，κ = 1.001 を用いて σ = (κ − 1) + max(diag( κ2 S T S +b HD
bh + Db HDb v ))) と，τ = (κ − 1) + 0.01 × ( κ + 2 ) ×
µ(D
h

v

2

κ

b †H S TS Db † )), max(diag(Db †HS TS Db † ))])
max([max(diag(Dvvhhのように設定した．ここで，diag は正方行列の対角成分を抜き出す操作である．実験結果を表 1 に示す．また，再構成画像を図 3 に，拡大図を図 4 に示す．従来法と比較して，再構成画像の PSNR が約 1.3 dB 向上し，計算時間も約半分に低減していることが分かる．実際に再構成結果
を見てみると，提案法によって再構成された MR 画像は，
濃淡の移り変わりがはっきりしていることが確認できる．
次に，T2 強調画像を参照画像 u として，T1 強調画像を間引かれた観測信号 y から再構成した．提案法の各パラメータは，変更せずに同じ設定を用いた．実験結果を表 2に示す．従来法と比較して，PSNR が約 0.8 dB 向上した．

T1 強調画像は解剖学的構造の描出に長けているため，T1
強調画像を参照画像として T2 強調画像を再構成する表 1の場合の方が，提案法がより有効性を発揮すると考える．

5

おわりに
本研究では，高速 MR 撮像で得られた間引き観測データ

から，エッジ情報と LiGME モデルを用いて MR 画像を高精度に再構成する手法を提案した．Ehrhardt らが提案した[1] S. W. Young, Nuclear Magnetic Resonance Imaging: Basic
Principles. New York: Raven Press, 1984.
[2] P. G. Morris, Nuclear Magnetic Resonance Imaging in Medicine and Biology. Oxford: Clarendon Press, 1986.
[3] D. J. Brenner and E. J. Hall, “Computed tomography—An
increasing source of radiation exposure,” New Engl. J. Med.,
vol. 357, no. 22, pp. 2277–2284, 2007.
[4] M. Lustig, D. L. Donoho, and J. M. Pauly, “Sparse MRI: The
application of compressed sensing for rapid MR imaging,”
Magnet. Reson. Med., vol. 58, no. 6, pp. 1182–1195, 2007.
[5] Y. Chen, X. Ye, and F. Huang, “A novel method and fast algorithm for MR image reconstruction with significantly undersampled data,” Inv. Prob. Imag., vol. 4, no. 2, pp. 223–240, 2010.
[6] S. Ravishankar and Y. Bresler, “MR image reconstruction fromhighly undersampled k-space data by dictionary learning,”
IEEE Trans. Med. Imag., vol. 30, no. 5, pp. 1028–1041, 2011.
[7] C. Tang, N. Inamuro, T. Ijiri, and A. Hirabayashi, “Compressedsensing MRI using double sparsity with additional training
images,” in Proc. ICASSP, New Orleans, 2017, pp. 801–805.
[8] M. J. Ehrhardt and M. M. Betcke, “Multicontrast MRI reconstruction with structure-guided total variation,” SIAM J. Imaging Sci., vol. 9, no. 3, pp. 1084–1106, 2016.
[9] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
[10] I. Selesnick, “Sparse regularization via convex analysis,” IEEE
Trans. Signal Process., vol. 65, no. 17, pp. 4481–4494, 2017.
[11] J. Abe, M. Yamagishi, and I. Yamada, “Linearly involved generalized Moreau enhanced models and their proximal splitting algorithm under overall convexity condition,” Inv. Prob.,
vol. 36, no. 3, 36 pages, 2020.
[12] I. Selesnick, A. Lanza, S. Morigi, and F. Sgallari, “Non-convex
total variation regularization for convex denoising of signals,”J. Math. Imaging Vis., vol. 62, no. 6–7, pp. 825–841, 2020.
[13] L. Condat, “A primal–dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms,” J. Optim. Theory Appl., vol. 158, no. 2, pp. 460–
479, 2013.
[14] L. Condat, D. Kitahara, A. Contreras, and A. Hirabayashi,“Proximal splitting algorithms: A tour of recent advances,with new twists,” preprint arXiv:1912.00137, 46 pages, 2019.

(a) 真の T2 強調画像(b) 逆離散フーリエ変換による初期画像 (22.81 dB)(c) 従来法による再構成結果 (26.78 dB)(d) 提案法による再構成結果 (28.09 dB)

図 3: 圧縮率 1/7 程度の観測信号からの T2 強調画像の再構成結果 (PSNR)(a) 真の T2 強調画像(b) 従来法による再構成結果(c) 提案法による再構成結果

図 4: 図 3 の再構成結果の拡大図