「MRIシミュレータを用いた独習パルスシーケンス〔標準編〕」でMRIをシミュレーション

今年はブロガーになる！と宣言したものの，1月の記事以降何も書いてなかった･･･．COVID-19のせいで幸か不幸かいろいろ仕事が無くなって時間ができたのでやっと更新できました．今回はMRIネタを持ってきました．先日，非常に面白い本を買ったのです．

元筑波大学の巨瀬先生が執筆された，MRIをPC上でシミュレーションするという新しい試みの本です．以前より，MRIシミュレータに関する研究をされておられ，その研究結果をこのような形で誰でも無料で使えるようにしていただいたようです．一応，磁気共鳴専門技術者である私は，このようなMRIに関するツールは大好きです．MRIの更なる理解，または後輩への教育等に使えれば良いなあと思い，このソフトを色々と試行錯誤して使ってみました．そこで今回はMRIをシミュレートするソフトBlochSolverの使用方法と雑感を報告したいと思います．

BlochSolverのダウンロード

まずは巨瀬先生の会社のホームページ上のソフトをダウンロードします．以下のURLからBlochSolver-CPUをダウンロードします．

http://mrisimulations.com/cpu_version_download/

同URLにマニュアルもありますので，これもダウンロードしておくべきでしょう．ダウンロード後，zipファイルを解凍しておきます．

スピンエコーをシミュレーション

ここから具体的にシミュレーションしていきたいと思います．まずは，コンベンショナルスピンエコーのTR・TEを変化させてシミュレーションしていきたいと思います．

1. TR・TEの設定

ダウンロードしたblochsolver-cpuフォルダ内に128_sequenceフォルダがあります．このフォルダの中にはマトリクスサイズ128×128のMRI画像のシーケンスシミュレーション用Pythonコードが書かれたファイルが入っています．今回はその中の2D_SpinEcho_128.seq.pyを編集してスピンエコーをシミュレーションします．
ちなみにこの.pyファイルはVScodeで編集することがマニュアル中で推奨されています．VScodeはエディタと言われるもので，プログラムを作るのに必要なものってイメージで大丈夫だと思います．VScodeは拡張機能を追加でき，非常に便利で，僕も研究でよく使っています．詳しいダウンロードの方法はググってください．
このファイルをVScodeで開いて5，6行目の箇所を以下の様に変更し保存します．

TR = 400.0e+3 # [us]
TE = 10.0e+3 # [us]

こうすることでTR=400ms，TE=10msのT1強調画像をシミュレーションできます．

2. command.txtの設定

TR・TEを設定した後，フォルダ内のcommand.txtファイルを編集します．この.txtファイルはwindowsであればメモ帳で開くことができます．VScodeでも開くことができます．このファイルは実際のシミュレーションをする際の様々な命令が書いてあります．以下のように書き直します．

set dim 128 128 128 
set subvoxel-dim 1 1 1
set actual-dim 256.0 256.0 256.0
set pd "./phantom/BrainPhantom_128/pd.flt"
set t1 "./phantom/BrainPhantom_128/t1.flt"
set t2 "./phantom/BrainPhantom_128/t2.flt"
set sequence "./128_sequence/2D_Spinecho_128.seq.py"
start
var PREFIX "./result/${DATE}/${SEQUENCE_NAME}-${SERIAL_ID}"
save complex "${PREFIX}-complex.flt"
quit

デフォルトから変更した箇所は1行目の
set dim 128 128 128
それと4～6行目の
set pd "./phantom/BrainPhantom_128/pd.flt"
set t1 "./phantom/BrainPhantom_128/t1.flt"
set t2 "./phantom/BrainPhantom_128/t2.flt"
最後に7行目の
set sequence "./128_sequence/2D_Spinecho_128.seq.py"
になります．1行目でマトリクスサイズを指定し，4～6行目で使用するファントムデータのプロトン密度値，T1値，T2値を指定しています．7行目で先ほど編集したシーケンスファイルを指定しています．

3. シミュレーションの実行

blochsolver.batをダブルクリックし実行します．以下の画像のような黒い画面が立ち上がり，処理が始まります． 処理が完了するとPress any key to continue . . .と表示されるのでなにかしらキーボードを押せば勝手にウィンドウが消え，処理が終わります．

4. 画像の確認

シミュレーションの結果をGUI_reconstruction_program//128//Reconstruction_and_display_for_128x128_2D_Cartesian.exeを使用して可視化します．このファイルを実行すると以下のウィンドウが立ち上がります．

左上のopenを選択し，resultフォルダ内に先ほどシミュレーションした結果のファイルが保存されているので，これを選択します．

シミュレーション画像ができました．一応専門家がこの画像を見た感想としては，画像の前後が逆になってる，側脳室付近に変なアーチファクトがでてる，T1強調画像というよりはプロトン密度強調画像っぽいといった感じです．一度画像を開くとresultフォルダ内に4つの画像が.fltファイルとしてそれぞれ別画像で保存されます．
なお作成された画像ファイルはImageJで開くことができるとマニュアルに書いてあり，実際にImageJではFile > Import > Raw...で以下の画像の設定で正しく読み込むことができました．

5. 別条件でシミュレーション

TR/TE=4000/100 msのT2強調画像も同様に作成し，ImageJで前後を変更してみました．

T2強調画像はコントラストも大丈夫そうでアーチファクトも見られません．T1強調画像のコントラストやアーチファクトについては，また本を読んで原因を解明したいと思います．

2020.7.17追記

このアーチファクトについて巨瀬先生からコメントを頂きました．このアーチファクトは短いTRが原因のflash band artifactのようです．しかし，臨床ではTR=400 ms程度でこのようなアーチファクトはでない印象があります．このファントムのT1値が近似により求められているため（詳細は本を参照），もしかしたら生体とこのファントムのT1値が異なっていることに一因があるのかもしれません．もしくは臨床機ではflash band対策をしたシークエンスがデフォルトになっているのかも？

まとめ

MRIは撮像する人により画質が大きく変化するモダリティであり，良い画像を取得する為には知識と経験が必要です．BlochSolverを利用すれば，実際にMRI装置を使用せずとも経験値を積むことができます．臨床で働く診療放射線技師にも非常に有用なソフトであると思います．今回はスピンエコーだけ試しましたが，他のシーケンスもシミュレーションできるようです．これはかなりおもしろいことができそうな予感･･･．

第7回岡山県西部医用画像研究会深層学習部会 –父と母どちらに似てる？–

はじめに

本記事は第7回岡山県西部医用画像研究会深層学習部会の参考資料となります．
これまで岡山県西部医用画像研究会深層学習部会では「ゼロから作るDeep Learning」を利用して深層学習の勉強をしました．前回の第6回では「ゼロから～」の第7章『畳み込みニューラルネットワーク(Convolutional Neural Network:CNN)』を勉強しました．約1年かけて当初の目標であったCNNの内容をとりあえず知ることができました（完璧に理解したとは言っていない･･･）．そこで今回，代表世話人である私杉本が現時点まで勉強したことを応用し，実際にプログラムを組み画像認識を行ってみました．せっかくここまで勉強したんだから，なにかアウトプットがしたかったのです．
ということで以下に第7回岡山県西部医用画像研究会深層学習部会として記事の内容を書いています．なるべく正確にしたつもりではありますが，初心者であるため間違いがあるかもしれません．ご容赦くださいm(__)m

目次

何をしたのか

子どもの顔が両親のどちらの顔に似ているかを客観的に判断できるAIを作成しました．ここでは私の子どもが私と妻，どちらのほうに似ているかを判断するAIを作りました．親となった私は，わが子を他人に見せた時に私に似てるといわれることが多いです．しかしこれがお世辞なのか，あるいは本当にそうなのか，それを客観的に知ることはできません．科学的にこの事実を判断したいという情熱からこのAIは作られました．またくだらないことをしてしまった．ほんとは医用画像を使って何かできればよかったのですが，なかなか面白そうな内容を思いつかなかったのでできませんでした．

その１：データ収集

第1回深層学習部会で少しお話したのですが，私の最近の趣味はライフハック．私は生活をより効率化することに喜びをおぼえる人です．私はGoogle photoを使ってスマホで撮った写真を自動的でバックアップしています．

photos.google.com

これがまたすごいんです．特にスマホの操作をせずに写真のバックアップをしてくれて，さらにそのバックアップした写真の顔を認識してその人ごとにアルバムを自動で作ってくれます．しかも無料．今回はその機能を利用して私と妻と子どものアルバムをダウンロードしました．

その２：顔のみを切り出す

ここからプログラムを書いていきます．実際にAIで人の顔の判断をさせるのに最も効率が良いのは顔のみの写真で学習させることです．今回はPython，OpenCVを使い顔の切り出しを行いました．以下，コードを書いています．

import sys, os
import glob
import cv2
import numpy as np
import matplotlib.pyplot as plt
import shutil
from PIL import Image

def filelist(path):
    return [p for p in glob.glob(path + '\\*')]
    
dir = '..\\rawdata'
file = filelist(dir)

data_path = '..\\data'
cropped_path = '..\\cropped'

for f in file:
    base, ext = os.path.splitext(f)
    if ext == '.jpg':
        shutil.copy(f, data_path)
    else:
        pass

#入力ファイルのパスを指定
in_jpg = data_path
out_jpg = cropped_path
 
#リストで結果を返す関数
def get_file(dir_path):
    filenames = os.listdir(dir_path)
    return filenames
 
pic = get_file(in_jpg)
 
for i in pic:
    # 画像の読み込み 
    image_gs = cv2.imread(in_jpg + '\\' +i)
 
    # 顔認識用特徴量ファイルを読み込む --- （カスケードファイルのパスを指定）
    cascade = cv2.CascadeClassifier("..\\opencv-master\\data\\haarcascades\\haarcascade_frontalface_alt.xml")
    
    # 顔認識の実行
    face_list = cascade.detectMultiScale(image_gs,scaleFactor=1.1,minNeighbors=1,minSize=(1,1))
    
    # 顔だけ切り出して保存
    no = 1;
    for rect in face_list:
        x = rect[0]
        y = rect[1]
        width = rect[2]
        height = rect[3]
        dst = image_gs[y:y + height, x:x + width]
        save_path = out_jpg + '/' + 'out_('  + str(i) +')' + str(no) + '.jpg'
        
        #認識結果の保存
        a = cv2.imwrite(save_path, dst)
        plt.show(plt.imshow(np.asarray(Image.open(save_path))))
        print(no)
        no += 1

上のプログラムが動くと顔だけ切り出されて保存されます．ここからが1番大変な作業でした．大量に切り出された顔をそれぞれフォルダ分けしていきました．もちろん手動．顔以外が認識されて保存されたりしていました．機械学習は前処理が8割と聞いたことがありましたが，まさにこのことかと思い知らされました．一週間くらいかかりました．なおデータは私の顔：785枚，妻の顔：421枚です．

その３：学習

私と妻の顔を切り出し，それらをCNNで学習させます．ちなみにGoogle colabを使っています．

colab.research.google.com

以下コードです．

# Google driveをマウント
from google.colab import drive 
drive.mount('/content/drive')

# ディレクトリの移動
%cd 'drive/My Drive/目的のフォルダ'

import os
import glob
import shutil
import numpy as np
import matplotlib.pyplot as plt

from tensorflow.python.keras import layers
from tensorflow.python.keras import models
from tensorflow.python.keras import optimizers
from tensorflow.python.keras.preprocessing.image import ImageDataGenerator


dir_train1 = './data/train/私'
dir_train2 = './data/train/妻'
dir_validation1 = './data/validation/私'
dir_validation2 = './data/validation/妻'

d_train1 = glob.glob(dir_train1+'/*')
d_train2 = glob.glob(dir_train2+'/*')
d_validation1 = glob.glob(dir_validation1+'/*')
d_validation2 = glob.glob(dir_validation2+'/*')

base_dir = './data'

train_dir = os.path.join(base_dir, 'train')
validation_dir = os.path.join(base_dir, 'validation')

train_datagen = ImageDataGenerator(rescale=1./255)
valid_datagen = ImageDataGenerator(rescale=1./255)

train_generator = train_datagen.flow_from_directory(
    train_dir,
    target_size=(256, 256),  # モデルの入力サイズ
    color_mode='rgb',  # 読み込み形式
    batch_size=20,  # バッチサイズ
    class_mode='binary')
valid_generator = valid_datagen.flow_from_directory(
    validation_dir,
    target_size=(256, 256),  # モデルの入力サイズ
    color_mode='rgb',  # 読み込み形式
    batch_size=20,  # バッチサイズ
    class_mode='binary')

# モデルの定義
np.random.seed(0)
model = models.Sequential()
model.add(layers.Conv2D(16, (3, 3), activation='relu', input_shape=(256, 256, 3)))
model.add(layers.MaxPooling2D((2, 2)))
model.add(layers.Conv2D(32, (3, 3), activation='relu'))
model.add(layers.MaxPooling2D((2, 2)))
model.add(layers.Flatten())
model.add(layers.Dense(256, activation='relu'))
model.add(layers.Dense(1, activation='sigmoid'))       

# モデルのコンパイル
model.compile(loss='binary_crossentropy',
              optimizer=optimizers.adam(),
              metrics=['acc'])

model.summary()

history = model.fit_generator(train_generator,
                              steps_per_epoch=100,
                              epochs=10,
                              validation_data=valid_generator,
                              validation_steps=50)

# モデルの保存
model.save('./models/model_01.h5')

今回作ったネットワークのモデルは以下になります． 畳み込み層プーリング層×2と全結合層×2のCNNになります．データを訓練用とテスト用に分割し，10epoch計算させました． 学習は以下のように行われました． これで学習ができました．訓練用データに対しては100%，テストデータに対しては93％の精度がでています． 以下のコードで可視化できます．

acc = history.history['acc']
val_acc = history.history['val_acc']
loss = history.history['loss']
val_loss = history.history['val_loss']

epochs = range(1, len(acc) + 1)

plt.plot(epochs, acc, 'bo', label='Training acc')
plt.plot(epochs, val_acc, 'b', label='Validation acc')
plt.title('Training and validation accuracy')
plt.legend()
plt.figure()

plt.plot(epochs, loss, 'bo', label='Training loss')
plt.plot(epochs, val_loss, 'b', label='Validation loss')
plt.title('Training and validation loss')
plt.legend()

plt.show()

その４：父と母どちらに似てる？

ここからが本番です．その３までで私と妻を判断するネットワークを作ることができました．このネットワークの中に子どもの写真を入力すると，私と妻のどちらかだとネットワークが判断します．

model = load_model('./models/model_01.h5')

img_predict = []

for image_name in os.listdir( './子どもの画像のフォルダ'):
    try:
        img = Image.open(os.path.join( './子どもの画像のフォルダ', image_name))
        img = img.convert("RGB")
        img = img.resize((256, 256))
        img_array = np.asarray(img)
        img_predict.append(img_array)
    except Exception as e:
        pass

img_predict = np.asarray(img_predict)

result_predict_classes = model.predict_classes(img_predict)

father = 0
mother =0

predict = result_predict.flatten()

for i in predict:
    if i ==0:
        mother += 1
    else :
        father += 1

print(father ,mother)

ここでは100枚の子どもの写真を入力しました．

結果

子どもの画像を入力した結果、
私：71，妻：29
に分類されました．つまり，93％の精度で私と妻を分類できるAIは，2人の子どもの顔を約70%の割合で私だと判断したという結果を導き出しました．やはり父親似でした．

まとめ

以上，くだらないモデルを作成することができました．1年間深層学習について勉強してきてとりあえずここまでたどり着くことができました．このモデルを作成する過程で本には書いていないわからないことがたくさん出てきて，深層学習の難しさをより味わうことができました．

最後に

プログラミングも深層学習もわからないことだらけですが，やらなければわからないことはわからないままです．昨今の人工知能ブームで医療の世界も変わりつつあります．さらにこれからはプログラミング教育も始まり，ますますこの分野から逃れることは難しくなると考えられます．診療放射線技師の業務に直接関わるところは少ないかもしれないですが，アプリケーションとして触れる機会は今後必ずやってくるでしょう．私はそれを使いこなせるように今後も深層学習をはじめ，機械学習についてボチボチ勉強し続けていきたいと思います．

ということで以上となります．明日からまた仕事頑張りましょう．