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

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

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

「対話形式だからわかりやすい!MRIの教科書」は初めてのMRIの教科書としてオススメ

本記事は「対話形式だからわかりやすい!MRIの教科書」の第1章と第2章の書評記事です。 本書は難しい数式を使わず噛み砕いた表現を使うことで、MRI の原理を気軽にわかりやすく理解することを目的に企画されたようです。一貫して対話形式で進む登場人物たちのやり取りを読むことで、難しいMRIの基礎原理を習得することができます。そのため本書はこれからMRIを学びたいと考えている人にとっての「最初の1冊」として適した本であると思います。

目次

きっかけ

ちょっと前にMRIの情報収集をしているときに面白そうなタイトルの本書を見つけました。タイトルを見て中身を読んでみたいと思いましたが、買おうとまでは思えませんでした。著者情報の記載はあるも、誰だかわからないし、個人出版のようなので情報に信頼性があるかがわからなかったからです。

そんなときに以下のツイートを見ました。Kindleストアで著者のTwitterアカウント(@MRI_taiwa)が記載されていたので、フォローしていたので気づくことができました。

このツイートに乗っかれば、気になってた本が読める!と思い、ツイートに返信してみました。数回のやり取りの後、この記事の執筆にいたりました。つまり、ブログでこちらの本を宣伝する代わりに本を読ませていただく(報酬を頂く)こととなりました。

可能な限りニュートラルな立場で書評をしたつもりですが、この点を理解した上で読んでいただければと思います。それでは以下に書評の内容を書いていきます。

本書はMRI基礎原理を学ぶために企画されている

この本の対象はMRIの基礎原理について興味がある、しかし、これまでじっくりと取り組んだことのない方、あるいは一度取り組もうとしたけれど挫折をした方だそうです。たしかに本書の内容はMRIの基礎について広く、浅く、易しく書かれています。

第1章はMRI検査全体についての大まかな説明

第1章を読むとMRIについての全体像を把握することができます。MRIの原理は非常に難しいため、まずは全体を知って、それから細かなところを深堀するほうがいいという筆者の考えから、この章を設けたようです。MRIのことをよく知らない人でも1章に書いてあることはなんとなく聞いたことがあると思います。

第2章では本格的にMRI装置のことについて

第2章ではMRI装置本体の原理を知ることができます。MRIは複数のコイルから成り立つという話から始まり、その後はそれぞれのコイルの役割について詳しく説明されています。

第2章は図をふんだんに使われており、筆者の熱量を感じられます。傾斜磁場の説明ではかなり頑張って作図されているということが垣間見えます。わかりやすく説明するために工夫を凝らしてある、素晴らしい本だと思いました。

本書の良い点

  • 対話形式は内容が頭に入ってきやすい
  • 動画(GIF)が面白い

本書はAさんとBさんと博士の3人の対話で内容が進んでいきます。僕が知る限り対話形式のMRIの本はありません。あえていうとすれば、荒木力先生の「MRI完全解説」をストーリーにしたものといった感じでしょうか。「MRI完全解説」はQ and A方式の教科書であり、趣向としては似ているのかなと思いました。ちなみに「MRI完全解説」は僕が磁気共鳴専門技術者試験に合格するために、最も参考にした書籍です。

また、本書は電子書籍の特徴を活かし、要所要所に動画(GIF)のリンクが貼られています。

これは非常に面白い試みだと思いました。MRIの原理は静止画だと理解が難しいため、動画で見ることができればその一助になると思います。これからもGIFをたくさん作っていただけるとさらに面白い本になると思いました。

本書の惜しい点

  • 登場人物がわからなくなる
  • 章ごとで販売されているため、1冊で考えたときの値段が高い

本書はAさんとBさんと博士の3人のやりとりで話が進んでいきます。僕は頭のメモリが弱いので、たびたびAさんとBさんがどっちがどっちかわからなくなりました。名前で表記していただければそのストレスもなかったのかなと思いました。

そして本書の値段ですが1章が300円、2章が500円となります。たしかにそれぞれの章で話は完結していますが、すべての章を合わせた1冊の値段で考えるとそれなりの値段になります。まとめ買いをすれば安く買える仕組みがあってもいいかと思いました。

しかしながら、本書はkindle unlimitedに登録すれば無料で読めます。30日間無料体験キャンペーンを長いことしているようです。本書が気になるけど、お金を出したくない方はkindle unlimitedも一つの方法かと思います。

今後の章に期待大!

2020年8月現在、まだ2章分しか販売されていませんが、その内容はMRIを学ぼうと考える人にとっての「初めの1冊」として非常に良いものだと思います。

とはいうものの、僕は本書に記載されていた「シールドフィンガー」や「シム」の言葉の由来について全く知りませんでした。臨床で直結する知識ではなく、いわゆる豆知識ですが、こういったセレンディピティともいえる発見があったので、結果的に僕は読んでよかったと思いました。

今後も続編発刊が予定されているようですので、続編にも大きく期待したいと思います。

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は2つ存在する

そもそも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

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

# 要素が0の入れ物を作っておく
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空間信号のみで画像再構成を行った場合、先に述べたように負の周波数成分があるため2つの画像が出現する。本例ではその一部が重なりあっている。重なりをなくすためにはk空間信号の最高周波数を下げる、すなわち空間分解能を下げる必要がある。これに対して複素k空間信号による再構成画像は負の周波数成分が発生しないため画像は1つのみである。したがって、複素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画像にどのような影響が出るかシミュレーションしてみたいと思います.

参考文献

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

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

決定版 MRI完全解説 第2版

決定版 MRI完全解説 第2版

  • 作者:荒木 力
  • 発売日: 2014/03/27
  • メディア: 単行本

Excelによる医用画像処理入門

Excelによる医用画像処理入門

被ばく線量情報抽出ソフト『DoNuTS』を作りました

f:id:radmodel:20200407234239p:plain

2020年4月1日より医療法施行規則の一部を改正する省令が施行され,CTなど一部の放射線科領域の医療機器を有する施設は被ばく線量の記録・管理が義務付けられました.多くの大規模施設では線量管理ソフトと呼ばれる高価なソフトを導入し,PACSと紐付けるなどして,簡便に線量の管理を行えるような体制を整えつつあります.しかし,大規模施設と比べて放射線科の検査件数が少ない中小規模施設では,新たに線量管理ソフトを導入することは困難であると考えられます.

そこでRDSR(Radiation Dose Structured Report)とPET画像のヘッダーからCTとPET/CTの被ばく線量情報を抽出し,csvで出力できるフリーソフトDoNuTS(Dose Counter for Nuclear Medicine and CT Systems)』を作りました.DoNuTSを使用することで,線量管理ソフトの導入が困難な施設においても,線量管理を簡便に行うことができると考えられます.今回はDoNuTSの紹介をしたいと思います.

なお,ソースコードは公開しており,無料でお使いいただけます.ただし,努力して作りましたがプログラミングに関しては素人が作ったため,バクがある可能性があります.ご使用は自己責任でお願いいたしますm(__)m

作成した方法

使用したプログラミング言語Python(version 3.6.9)です.windows上でソフトとして使用する為には.pyファイルを.exeファイルにする必要があります.今回はpyinstallerを使用してexe化しました.

ダウンロードの方法

ソフトは以下のリンクからダウンロードできます.

リンクをクリックすると,僕のGithubのページに飛びます.

f:id:radmodel:20200407224432p:plain

右上のあたりにある緑色のClone or downloadをクリックすると,ファイルのダウンロードが始まります.ダウンロードしたファイルの中にsrcという名前のフォルダがあります.そのsrcの中にdistというフォルダがあり,その中にDoNuTS.exeというファイルがあります.それがDoNuTSの本体となります.

使用方法

データの準備

まずはRDSRとPET画像をDoNuTSを動作させるPC上にexportしておく必要があります.このとき,1つのフォルダ内にデータを取得したいRDSRおよびPET画像を保存してください.なお,DoNuTSはRDSRとPET画像のみを探して処理を行うため,他のファイルがフォルダ内に存在していても動作に問題はありません.

DoNuTSの実行

データの準備ができたらDoNuTSをダブルクリックして実行します.するとコマンドプロンプトが立ち上がります.黒い画面のアレです.次に,取得したいデータが入ってあるフォルダを指定するよう指示されるので,指示に従い先ほど準備したフォルダを選択します.処理が終わると以下の図のようなポップアップが表示されます.

f:id:radmodel:20200407230623p:plain

これで処理は完了です.なお,出力されるデータは

  • Mean CTDIvol
  • DLP
  • X-Ray Modulation Type
  • CTDIw Phantom Type
  • Acquisition Protocol
  • Target Region
  • CT Acquisition Type
  • Procedure Context
  • Exposure Time
  • Scanning Length
  • Exposed Range
  • Nominal Single Collimation Width
  • Nominal Total Collimation Width
  • Pitch Factor
  • Identification of the X-Ray Source
  • KVP
  • Maximum X-Ray Tube Current
  • Mean X-Ray Tube Current
  • Exposure Time per Rotation
  • Device Manufacturer
  • Device Serial Number
  • Manufacturer Model Name
  • Patient ID
  • Study Date
  • Patient Name
  • Study Description
  • Patient BirthDate
  • Patient Sex
  • Patient Age
  • Patient Size
  • Patient Weight
  • CT Dose Length Product Total
  • Radionuclide Total Dose (PET画像が存在した場合)

2020.7.17 追記

(出力できるデータが増えたため,リストを変更しました.)

以上です.出力されたcsvは,先ほど指定したフォルダ内に保存されています.

実際の業務における使用例

診療放射線技師会や医学放射線学会のホームーページで線量管理に関するガイドラインが発表されています.その中身を確認してみるとおおよそ以下のことをしておけば問題はないと考えられます.

  1. 線量の記録
    • dose summaryをPACSに転送しておく
  2. 線量の管理
    • 記録してある線量情報をサンプリングし,自施設の線量を算出
    • DRL等を利用して自施設の線量を評価し,最適化を行う

ここでDoNuTSは線量の管理で効果を発揮します.たとえばCTの線量管理の場合,検査終了後に撮影した患者のRDSRをPACS等に保存する運用をしていき,線量の管理を行う際に保存しておいたRDSRをexportし,DoNuTSを使用します.そうすることでDLPや体重が患者ごとにcsvとして出力されるので,Excel等をうまく利用すれば被ばく線量の最適化を簡便に行うことができると考えられます.

まとめ

最近の装置はRDSRを出力できるのが一般的になっています.RDSRは構造化されたレポートという名前の通り,データを取り出すことさえできれば,非常に便利なファイルです.しかし,実際は高価な線量管理ソフトを導入しなければ,RDSRからデータを取り出すこと自体困難です.DoNuTSはRDSRの中の線量管理に必要な情報を簡便に取り出すことができるので有用だと思います.是非使ってみて下さい.

また,ソフトに関する要望やバグを発見した場合は,Github上のissueに書き込むか,メール(deepokayama@gmail.com)をしていただくか,なんらかの手を使って僕にお知らせいただければ幸いです.

「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をダブルクリックし実行します.以下の画像のような黒い画面が立ち上がり,処理が始まります. f:id:radmodel:20200310011908p:plain 処理が完了するとPress any key to continue . . .と表示されるのでなにかしらキーボードを押せば勝手にウィンドウが消え,処理が終わります.

4. 画像の確認

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

f:id:radmodel:20200310012730p:plain

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

f:id:radmodel:20200310013141p:plain

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

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

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

f:id:radmodel:20200310015255p:plain

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:データ収集

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

photos.google.com

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

その2:顔のみを切り出す

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

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枚です.

その3:学習

私と妻の顔を切り出し,それらを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')

今回作ったネットワークのモデルは以下になります. f:id:radmodel:20200105140248j:plain 畳み込み層プーリング層×2と全結合層×2のCNNになります.データを訓練用とテスト用に分割し,10epoch計算させました. 学習は以下のように行われました. f:id:radmodel:20200105140219j:plainこれで学習ができました.訓練用データに対しては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()

f:id:radmodel:20200105140625p:plainf:id:radmodel:20200105140637p:plain

その4:父と母どちらに似てる?

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

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年間深層学習について勉強してきてとりあえずここまでたどり着くことができました.このモデルを作成する過程で本には書いていないわからないことがたくさん出てきて,深層学習の難しさをより味わうことができました.

最後に

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

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