放射線技術×データサイエンス

放射線技術×データサイエンス

診療放射線技師がPythonとRとImageJを使って色々します

Motion Artifact Generator (MAGE) をPythonで実装した

楽しみにしていた以下の論文を読みました。

www.jstage.jst.go.jp

この論文はMAGE (motion artifact generator)と呼ばれるMRI画像のモーションアーチファクトジェネレータを開発した論文です。アーチファクトを擬似的に作成できれば、元の画像を教師画像とした深層学習モデルを作成することができます。MAGEはモデル作成時に必要なモーションアーチファクト画像を大量に作成するためのシステムだそうです。

MAGEは以下の図のようにしてモーションアーチファクト画像を作成します。MRI画像を複数のパターンでシフトさせ、それらの画像をフーリエ変換して得られた擬似k-spaceを組み合わせることでモーションアーチファクトを生成しているようです*1。詳細は論文をご覧ください。

塚本、室(2021)から図を引用

面白そうだと思ったので、MAGEをpythonで実装してみました。今回のブログではそれを報告したいと思います。コードはgithubにあげています。

github.com

目次

動作環境

動作の検証

使用画像

今回使用した画像は、エムアールアイシミュレーションズ様が開発したBlochSolverで以前シミュレーションした人体頭部のSEのT2WI画像です。

www.blochsolver.org

radmodel.hatenablog.com

まず、シミュレーションで得たデータを読み込みます。

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

# 作成したmageモジュールをインポート
import mage

# BlochSolverでシミュレートして作成したSE画像を読み込むための関数
def sim2img(path, matrix):
    with open(path, 'rb') as f:
        d = f.read()

    row_data = Image.frombytes('F', (matrix[0]*2, matrix[1]), d)
    row_data = np.array(row_data)
    row_data.shape
    img = np.zeros(matrix, dtype=np.complex64)
    for col in range(matrix[1]):
        for row in range(matrix[0]):
            img[row, col] = row_data[row, col*2] + (row_data[row, (col*2+1)]*1j)
    return img            

# BlochSolverで作成したMRI画像を読み込む
path = '2D_SpinEcho-3-complex.flt'
kspc = sim2img(path, [256,256])

# k-spaceを表示
mage.plot_kspace(kspc)

# 強度画像を表示
t2wi = np.abs(np.fft.ifftshift(np.fft.ifft2(kspc)))
plt.imshow(t2wi, cmap='gray')

MRI画像をシフトさせる

次にMRI画像をシフトさせてみます。今回は論文と同様に前後方向に-10〜10 pixel移動した画像を作成します。

# Image形式に変換
data = Image.fromarray(t2wi)

# SEでモーションアーチファクトを生成
# mageモジュールのMAGESEクラスからSE用MAGEのインスタンスを作成
mage_img = mage.MAGESE(data)

# シフトさせる最大のピクセル量
move_pixel = 10
move_range = np.arange(move_pixel*-1, move_pixel+1, 1)

# 前後方向に動かすメソッドmove_vertical()を実行
mage_img.move_vertical(move_pixel)

# 前後方向に動いた画像を表示
for k, _ in enumerate(move_range):    
    plt.subplot(5, 5, k+1)
    plt.axis('off')
    plt.imshow(mage_img.vertical_moved_image[k], cmap='gray')
    plt.gca().title.set_text('Vertical move ('+str(k-10)+' pixel)')

擬似k-spaceの作成

次にシフトしたMRI画像をフーリエ変換することで擬似k-spaceに変換し、それらの複数の擬似k-spaceを位相エンコード方向にランダムサンプリングして、一枚のk-spaceを作成します。

私の実装では先ほどMRI画像をシフトさせたときに同時にk-spaceも作成し、またシフトした擬似k-spaceデータからランダムサンプリングして1枚のk-spaceを作成しています。そのため、あえてk-spaceを作成するコードを書く必要はありません。

モーションアーチファクトMRI画像の作成

シフトした画像をランダムにサンプリングしてモーションアーチファクトを模擬した擬似k-spaceを逆フーリエ変換します。これでモーションアーチファクトを生じさせたMRI画像が作成できます。

# 強度画像にして表示
img = np.abs(mage.kspace2img(mage_img.concate_vertical_moved_kspace))
plt.imshow(img, cmap='gray')

以上です。微妙にモーションアーチファクトが出ていると思います。

モーションアーチファクトの強度は調整できます。MAGEのコードの以下の部分のsdを大きくするほどアーチファクトが強く出るようになります。

def create_random_kspace(self, move_range, moved_kspace):
    concat_moved_kspace = np.zeros(self.matrix, dtype=np.complex64)

    # 乱数を設定
    mean = np.median(move_range)
    sd = 5
    random_number = self._get_trunked_normal(mean, sd, 0, np.max(move_range), self.matrix[self.pe])
    random_number = np.ceil(random_number).astype(int)

TSEおけるMAGE

MAGEはSEとTSEの2種類のシークエンスでモーションアーチファクトを生成することができます。そこで私のpythonでの実装もMAGESEMAGETSEの2つのclassを作成しました。どちらのシークエンスでもアーチファクト画像を作成できます。さきほどはSE版でアーチファクト画像を生成しました。TSE版では同じシフト量で以下のような画像が作成できました。

# TSE用のMAGEのインスタンスを作成
# 2つ目の引数はecho train length
tse_img = mage.MAGETSE(data, 8)

# シフトさせる最大のピクセル量
move_pixel = 10

# 前後方向に動かすメソッドmove_vertical()を実行
tse_img.move_vertical(move_pixel)

# 強度画像にして表示
tse_mag = np.abs(mage.kspace2img(tse_img.concate_vertical_moved_kspace))
plt.imshow(tse_mag, cmap='gray')

論文の画像とはちょっと違うような気がしますが、なんとなくMAGEを再現できました。

*1:k-space上でモーションアーチファクトを生成する方法はtwitterのフォロワーさんから教えていただいたこちらの論文が参考になります(s0pさんありがとうございます)