Pythonを使ってk-spaceを理解する[フーリエ変換編] - 放射線技術×データサイエンス

k-spaceはよくこんな感じで説明されてます．

f:id:radmodel:20200510004035p:plain

k-spaceはMRIで受信した信号を格納したものと私は理解してます．k-spaceを逆フーリエ変換するとMRI画像が取得できます．
ところで，大手医療機器メーカーが販売しているMRI装置では，ユーザーは基本的にk-spaceを見ることはできません．生データが格納されているので，ユーザーに知られたくない情報が明らかになってしまう恐れがあるからだと私は思ってます（各社メーカー独自のフィルタ等の処理がされている？）．そのため，我々診療放射線技師はk-spaceの存在を知ってはいるものの，実際に見る機会は少ないです．k-spaceはアーチファクトやSNRや分解能を理解するのに不可欠なものです．もしそれを直接触ることができたら，MRIについてもっと理解が深まると思いますが，簡単にはできませんでした．

しかし，今はできるようになりました．前回紹介した巨瀬先生のMRIシミュレータ本ではk-spaceを取得することができます．

radmodel.hatenablog.com

前回はBlochSolverの簡単な使い方を試したのですが，今回はMRIシミュレータで作成したk-spaceをPythonを使ってフーリエ変換を行い，MRI画像を再構成してみました．装置がやることを自分で実際にやってみて理解を深めることが目的です．

k-spaceは２つ存在する

そもそもMRIは信号を直角位相感受性検波という手法を用いて取得しています．基準周波数と90度位相のずれた基準周波数の両方の信号を受信する方法です．この方法を用いることで正確に信号を認識できるらしいです．また，基準周波数の信号を実信号，90度ずれた信号を虚信号と言います．つまり，2種類の信号のデータが存在します．ということはk-spaceは実信号と虚信号の2つ分存在すると考えることができます．研究会や学会で演者の人が見せるスライドには，この記事の最初に示したような1つのk-spaceで説明されることが多い印象ですが，正確には2つデータがあるので実信号のk-spaceと虚信号のk-spaceの2つを考える必要があります．

実際にk-spaceを見てみる

ここからMRIをシミュレーションしていきます．シーケンスはスピンエコーを使用しました．シミュレーションの対象はスライス評価ファントムです．以下その画像．

f:id:radmodel:20200510012505p:plain

Matrix = 256×256，TR/TE = 4000/100 msecでシミュレーションしました．シミュレーションのやり方は前回の記事を参考にしてください．

さて，シミュレーションで作成されたデータはResultフォルタ内にシーケンス名-complex.fltという形で保存されています．このデータは，MRI信号の実信号と虚信号が32bit浮動小数点で交互に保存されているようです．まずはPythonを使って交互に保存されたデータをそのまま見てみます．

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

# データが存在するパスを指定
data = 'シーケンス名-complex.flt'

# データを読み込む
with open(data, 'rb') as f:
    d = f.read()

# バイナリデータをイメージ形式に変換
raw = Image.frombytes('F', (512,256), d)

# 画像を保存
plt.imshow(raw)
plt.savefig("raw.png")

f:id:radmodel:20200510015452p:plain

次に，このデータを実信号と虚信号に分けてみます．

# 要素が０の入れ物を作っておく
real = np.zeros([256,256])
imag = np.zeros([256,256])

# データを分割
for i in range(256):
    for j in range(256):
        real[i,j] = raw[i,2*j]
        imag[i,j] = raw[i,2*j+1]

# それぞれを保存
fig = plt.figure()

ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)

ax1.imshow(real, cmap='gray')
ax1.set_title('Real')
ax2.imshow(imag, cmap='gray')
ax2.set_title('Imaginary')

plt.savefig('real_imag.png')

f:id:radmodel:20200510020404p:plain

フーリエ変換

ここから診療放射線技師養成校で習うものの中でも鬼門となるフーリエ変換のパートに入ります．実信号と虚信号のデータを取得できたので，これらを使って絶対値表現されたMRI画像を作ります．

ここで私は今までたいへんな勘違いをしていることに気づきました．通常用いるいわゆるMRI画像は，実画像と虚画像との絶対値画像であるとは知っていましたが，実画像は実信号が入ったk-space，虚画像は虚信号が入ったk-spaceをそれぞれフーリエ変換して，それらの絶対値をとって作成するものだと思ってましたが，違いました．実際に実信号だけで実画像を作成してみます．

# ダメな例
bad_real = np.fft.ifft2(real)
bad_real = np.fft.ifftshift(bad_real)
bad_real = bad_real.real
plt.imshow(bad_real, cmap='gray')
plt.savefig('bad_real.png')

f:id:radmodel:20200510215623p:plain

実際にやってみてびっくり．実信号のみでフーリエ変換を行うとひっくり返った画像が重なってしまいます．
このことについて，金沢工業大学の樋口先生のホームページにわかりやすい説明がありました．

以下樋口先生のホームページからの引用です．

実k空間信号のみで画像再構成を行った場合、先に述べたように負の周波数成分があるため２つの画像が出現する。本例ではその一部が重なりあっている。重なりをなくすためにはk空間信号の最高周波数を下げる、すなわち空間分解能を下げる必要がある。これに対して複素k空間信号による再構成画像は負の周波数成分が発生しないため画像は１つのみである。したがって、複素k空間信号の方がフーリエ空間上の点を最大限利用できる。

こんなの知らなかった…でも，シミュレーションでk-spaceを触れたからこそ，新たに知ることができました．

話を戻して，ではどのようにすれば実画像と虚画像が作成できるかというと，樋口先生が記載しているように実信号と虚信号を複素数にして1つにまとめ，それをフーリエ変換すればいいようです．つまり，k-spaceを以下のように複素数で表現します．

$F(u)=F_{Real}(u)-\mathrm{i}F_{Imag}(u)$

ここで $F_{Real}(u)$ が実信号のk-spaceのデータ， $F_{Imag}(u)$ が虚信号のk-spaceのデータです．次に簡単のため，1次元の逆フーリエ変換の式を以下に示します．

$f(x)=\int^{\infty}_{-\infty} F(u)e^{i2\pi ux} du$

またオイラーの公式は以下の通り

$e^{\mathrm{i}2\pi ux}=\cos(2\pi ux)+\mathrm{i}\sin(2\pi ux)$

よって最初の逆フーリエ変換の式は以下のように変形できます．

$f(x)=\int^{\infty}_{-\infty} F(u)\cos(2\pi ux)du+\mathrm{i} \int^{\infty}_{-\infty} F(u)\sin(2\pi ux)du$

このとき，実画像と虚画像は以下のように表されます．

$f_{Real}(x)=\int^{\infty}_{-\infty} F(u)\cos(2\pi ux)du$

$f_{Imag}(x)=\int^{\infty}_{-\infty} F(u)\sin(2\pi ux)du$

これらを書き直すと，

$f(x)=f_{Real}(x)+\mathrm{i}f_{Imag}(x)$

となります．ここで $f_{Real}(x)$ が実画像， $f_{Imag}(x)$ が虚画像になります．そして通常よく用いられる絶対値MRI画像は以下の式で求められます．

$|f(x)|=\sqrt{f_{real}^2(x)+f_{Imag}^2(x)}$

実信号のk-spaceと虚信号のk-spaceの2つを複素数として結合して，逆フーリエ変換をすれば，実画像，虚画像が複素数として表現された正しい値が取得できるようです．この式変形は下に書いてある参考文献の『Excelによる医用画像処理入門』を参考にしました．どこか間違ってたら教えていただければありがたいです．

ということでこれをPythonでやってみます．

# 複素数の配列を作って実信号と虚信号を1つの複素数とする
com = np.zeros([256,256], dtype=np.complex64)
for i in range(256):
    for j in range(256):
        com[i,j] = complex(real[i,j],imag[i,j])

# 複素数データを逆フーリエ変換をする
ifft = np.fft.ifft2(com)
ifft = np.fft.ifftshift(ifft)

# 実画像，虚画像，絶対値画像に分ける
real_ifft = ifft.real
imag_ifft = ifft.imag
abs_ifft = np.abs(ifft)

# データを表示
fig = plt.figure(figsize=(8.0, 4.0))

ax1 = fig.add_subplot(1,3,1)
ax2 = fig.add_subplot(1,3,2)
ax3 = fig.add_subplot(1,3,3)

ax1.imshow(real_ifft, cmap='gray',vmin=0,vmax=1)
ax1.set_title('Real')
ax2.imshow(imag_ifft, cmap='gray',vmin=0,vmax=1)
ax2.set_title('Imaginary')
ax3.imshow(abs_ifft, cmap='gray')
ax3.set_title('Magnitude')

plt.savefig('final_data.png', )

f:id:radmodel:20200510032938p:plain

こうするとなぜか実画像と虚画像に縞模様が･･･．別々で保存してみます．

# real
plt.imshow(real_ifft, cmap='gray',vmin=0,vmax=1)
plt.title('Real')
plt.savefig('real_ifft.png')

# imag
plt.imshow(imag_ifft, cmap='gray',vmin=0,vmax=1)
plt.title('Imaginary')
plt.savefig('imag_ifft.png')

f:id:radmodel:20200510033140p:plain

f:id:radmodel:20200510033205p:plain

縞模様が消えました．なぜ縞模様がでるのだ…とりあえず絶対値画像が作れたのでよしとします．

まとめ

k-spaceってほんと難しい．直接見れないのもあるし，フーリエ変換とかいうよくわからない強敵もいるし，なにかと理解を拒む壁があります．でもシミュレーションデータで仮想的にk-spaceを作って自分で再構成をすることで，少し理解できた気になれました．まだまだ知らないことだらけです．次はもう少しフーリエ変換について勉強して記事を書こうと思います．そのあとはk-spaceをいろいろいじってMR画像にどのような影響が出るかシミュレーションしてみたいと思います．