# 京都大学がんプロセミナー(2019/2/17)
# 「Deep learningによる医用画像解析の深化と進化」
# オンライン演習用資料
# - Google Colaboratoryを利用した臓器セグメンテーション演習 -
#
# Copyright @2019 e-Growth Co., Ltd.
# http://www.egrowth.co.jp/
#
# このレッスンではDICOMスライス方向(アキシャル方向)ではなく、コロナル(正面断面方向)で学習してみることを試します
# データを既にダウンロード済みの場合はダウンロードの部分をスキップして問題ないです
# データのダウンロード
#
# LIDC-IDRIのCTデータに肺の輪郭を付加したものDICOM-RTデータに対し、
#「Growth RTV」でCT 3D Raw(入力)および肺領域の3D Raw(教師データ)を抽出したもの
!wget https://www.dropbox.com/s/yqcr2laea475rhg/data_set_3d.zip
# データの解凍
# 解凍後は data_set_3d という名のフォルダが生成される
!unzip data_set_3d.zip
!ls
### 学習用クラスのダウンロード ###
# スムーズな実習を進めるため、セグメンテーション用のクラスは予め実装しておきました
# 今回の実習はセグメンテーションの処理にU-Netと呼ばれるモデルを利用します
# U-Netの論文URLはこちら https://arxiv.org/pdf/1505.04597.pdf
!wget https://www.dropbox.com/s/rvnwmd6nbx6kvsj/seminor_inc.zip
!unzip seminor_inc.zip
!ls
import os, sys, glob
import numpy as np
import cv2
import keras
from keras.optimizers import Adam
import keras.backend as K
sys.path.append("seminor_inc/") #seminor_incフォルダ以下のソースコードも参照できるように設定
from my_functions import *
from unet import *
####### データの前処理 #######
dataDir = "data_set_3d/" #データフォルダ
#CTファイル一覧を取得
ctRawList = glob.glob(dataDir + "*ct.raw")
ctRawList.sort()
### データ読み込み ###
dataX = []
dataY = []
for rawPath in ctRawList:
infoTxt = rawPath.replace("ct.raw", "info.txt") #解像度を記録しているTxtファイル名に置き換え
shape = np.loadtxt(infoTxt,delimiter=",", dtype=np.int32) #解像度をint型として読み取る
shape = (shape[0], shape[1], shape[2], 1) #学習用画像及び教師データはともにグレースケールの1チャンネルであるので、チャンネル数情報として追加しておく
rawCT = readRaw(rawPath, dataType=np.int16) #CTの3D Rawを読み込む(この時点では１次元配列)
rawCT = rawCT.reshape( tuple(shape) ) #4次元配列(Z, Y, X, Channel)に戻す
print(rawPath, rawCT.shape) #rawCTのshapeを確かに4次元配列になっているかを確認してみる
rawLungPath = rawPath.replace("ct.raw", "Lung.raw")
rawLung = readRaw(rawLungPath, dataType=np.uint8) #同様にしてLung領域を示す3D Rawを読み込む. ただしCTと違い、uint8型(0～255)であることを指定
rawLung = rawLung.reshape( tuple(shape) ) #CTと同じ4次元配列に整形
rawCT = np.clip(rawCT, -1000, 1000) #CT値は範囲が広すぎるので、対象臓器が含まれる範囲にクリップ
dataX.append(rawCT) #CTを学習時の入力用配列に登録
dataY.append(rawLung) #CTを学習時の教師用配列に登録
if False:
### 試しに中央のスライスとそのLung領域を表示してみる
slicePos = int(rawCT.shape[1] / 2) #中間の位置を算出
sliceCT = rawCT[:, slicePos, :] #CTスライスを取り出す
sliceUint8 = getSliceByLut(sliceCT, wl=0, ww=2000) #WLとWWを指定して表示用画像に変換
sliceLung = rawLung[:, slicePos, :] #同じ位置のLung領域を取り出す
sliceUint8 = cv2.resize(sliceUint8.astype("uint8"), (256, 256))
sliceLung = cv2.resize(sliceLung.astype("uint8"), (256, 256))
img = np.hstack((sliceUint8.reshape((256, 256)), sliceLung.reshape((256, 256)))).astype("uint8")
import matplotlib.pyplot as plt
plt.imshow(img)
plt.gray()
plt.show()
exit()
### 学習用にCT値を -1から1の範囲に正規化
### これは好きな範囲で良い(例えば0から1の範囲に正規化しても良い)
dataX = np.array(dataX) / 1000.0
#dataX = np.clip(dataX, -1000, 1000) / 1000.0
### 教師用データを0～1の範囲に正規化
### 重要な注意点として、学習モデルの出力層の活性化関数に合わせて正規化する必要がある
### 例えば、sigmoid関数なら0から1、tanh関数なら-1から1の範囲に合わせる
### 今回はsigmoidを出力に使っているので0から1
dataY = np.array(dataY) / 255.0
### 学習用と検証用にデータを分割 ###
trainNum = int( len(dataX) * 0.8 ) #学習用データ数。今回は8割(10データ中8データ)を学習に使う
testNum = len(dataX) - trainNum #検証用データ数
if True:
trainX = dataX[:trainNum] # 学習用入力
trainY = dataY[:trainNum] # 学習用教師データ
testX = dataX[trainNum:] # テスト用入力
testY = dataY[trainNum:] # テスト用教師データ
else:
trainX = dataX[testNum:] #学習用入力
trainY = dataY[testNum:] #学習用教師データ
testX = dataX[:testNum] #テスト用入力
testY = dataY[:testNum] #テスト用教師データ
#各データのShapeを確認
print(trainX.shape)
print(trainY.shape)
print(testX.shape)
print(testY.shape)
#本当に正規化が正しく行われているかを最初のデータで見てみる
print("dataX[0]", np.min(trainX[0]), np.max(trainX[0]))
print("dataY[0]", np.min(trainY[0]), np.max(trainY[0]))
####### 学習の実施 #######
####### lesson2ではアキシャル方向で学習をためしたが、このlessonではコロナル(正面断面)方向で学習を試してみます ######
#モデルを読み込み
model = getUnet(256, 256, 1)
# ダイス係数を計算する関数
def dice_coef(y_true, y_pred):
y_true = K.flatten(y_true)
y_pred = K.flatten(y_pred)
intersection = K.sum(y_true * y_pred)
if( 0 == K.sum(y_true) + K.sum(y_pred) ):
return 1.0
else:
return (2.0 * intersection) / (K.sum(y_true) + K.sum(y_pred))
# 損失関数
def dice_coef_loss(y_true, y_pred):
return 1.0 - dice_coef(y_true, y_pred)
#損失関数はダイス係数(1.0より小さいほど損失大)、オプティマイザーはAdamを利用してコンパイル
model.compile(loss=dice_coef_loss, optimizer=Adam(lr=0.0005), metrics=[dice_coef])
# 学習したモデルの保存用フォルダを生成
modelDir = "model/"
if not os.path.isdir(modelDir):
os.makedirs(modelDir)
fromEpoch = 0 #開始エポック
totalEpoch = 3 #終了エポック
steps = 200 #1エポック内の学習回数
batchSize = 9 #1学習あたりのバッチサイズ
#学習済み重みをリロードして、続きから学習する場合
if False:
model.load_weights("model/ep_%08d_w.h5" % (fromEpoch-1))
#入力画像を(256,256)ピクセルに線形補間し、(256,256,1)に整形して出力する関数
def getResizeImage(img):
imgResize = cv2.resize(img, (256, 256)).reshape((256, 256, 1))
return imgResize
# totalEpoch - fromEpoch回のエポックを繰り返す
for iEpoch in range(fromEpoch, totalEpoch):
# エポックごとにsteps回の学習を実施
# つまり、steps x (totalEpoch - fromEpoch) 回の学習と重み修正が行われる
for iStep in range(steps):
trainBatchX, trainBatchY = [], []
dataIds = np.random.randint(0, trainNum, 3) #ステップごとに3つデータをランダムに決定
### データごとにスライスを任意に取得する処理
for dataId in dataIds:
#sliceNum = trainX[dataId].shape[0] #そのデータ内のスライス数を参照
#コロナル(冠状断面、正面断面)方向にデータをスライスしてみる
sliceNum = trainX[dataId].shape[1]
sliceIds = np.random.randint(0, sliceNum, 3) #当該データ内で3スライスをランダムに決定
for sliceId in sliceIds:
xTmp = getResizeImage(trainX[dataId][:, sliceId, :]) #コロナル方向ではデータごとに画像サイズが異なるので、これを256x256に線形補間
yTmp = getResizeImage(trainY[dataId][:, sliceId, :]) #コロナル方向ではデータごとに画像サイズが異なるので、これを256x256に線形補間
trainBatchX.append( xTmp ) #バッチに正規化済みスライスCTを入力として加える
trainBatchY.append( yTmp ) #バッチに正規化済みLung領域を教師データとして加える
# モデルに渡す際は必ずnumpy型に変換
trainBatchX = np.array(trainBatchX)
trainBatchY = np.array(trainBatchY)
print(trainBatchX.shape, trainBatchY.shape)
### 取得したバッチで学習実施
loss = model.train_on_batch([trainBatchX],[trainBatchY])
print("epoch", iEpoch, "step", iStep, "loss:", loss[0], "acc:", loss[1])
#学習した重みを保存
weightPath = "%s/ep_%08d_w.h5" % (modelDir, iEpoch)
model.save_weights(weightPath)
#エポックごとに学習が終わったら, テスト用データで検証してみる
testBatchX, testBatchY = [], []
### testデータごとにスライスを任意に取得する処理
for dataId in range(testNum):
sliceNum = testX[dataId].shape[0] # そのデータ内のスライス数を参照
#sliceIds = np.random.randint(0, sliceNum, ) # 当該データ内で3スライスをランダムに決定
for sliceId in range(sliceNum):
xTmp = getResizeImage(testX[dataId][:, sliceId, :]) #コロナル方向ではデータごとに画像サイズが異なるので、これを256x256に線形補間
yTmp = getResizeImage(testY[dataId][:, sliceId, :]) #コロナル方向ではデータごとに画像サイズが異なるので、これを256x256に線形補間
testBatchX.append( xTmp ) # バッチに正規化済みスライスCTを入力として加える
testBatchY.append( yTmp ) # バッチに正規化済みLung領域を評価用教師データとして加える
# numpy型に変換
testBatchX = np.array(testBatchX)
testBatchY = np.array(testBatchY)
# 評価実施
score = model.evaluate([testBatchX], [testBatchY])
print("Evaluate result:", iEpoch, "loss:", score[0], "acc:", score[1])
# 上記設定では僅か3エポックしか学習していないため、テストデータに対する抽出精度もまだ低くいので、
# 続きから追加で学習させてどこまで精度が上がるかを確認!
# 上記コードに対して以下のように変更して、再度学習をしてみる
# ■変更箇所 1
# fromEpoch = 0
# ↓
# fromEpoch = 3 # 3エポック目から追加学習
#
# ■変更箇所 2
# totalEpoch = 3
# ↓
# totalEpoch = 10 # 全部で10エポック分学習実施
#
# ■変更箇所3
# if False:
# model.load_weights("model/ep_%08d_w.h5" % (fromEpoch-1))
# ↓
# if True: # False から Trueに変更し、前回学習した重みを読み込む
# model.load_weights("model/ep_%08d_w.h5" % (fromEpoch-1))
#
# 以上で変更完了です。再度再生ボタンをおして、実行してみよう！
# 検証用データに対する予測結果
import matplotlib.pyplot as plt
model = getUnet(256, 256, 1)
#指定した重みを読み込み
wEpoch = 2
model.load_weights("model/ep_%08d_w.h5" % (wEpoch))
for iData in range(testX.shape[0]):
for slicePos in range(testX[iData].shape[1]):
sliceCT = testX[iData][:, slicePos, :] # CTスライスを取り出す
sliceCT = getResizeImage(sliceCT)
sliceLungPred = model.predict(sliceCT.reshape((1, 256, 256, 1))) * 255
sliceCT *= 1000
sliceUint8 = getSliceByLut(sliceCT, wl=0, ww=2000) # WLとWWを指定して表示用画像に変換
sliceLung = testY[iData][:, slicePos, :] * 255 # 同じ位置のLung領域を取り出す
sliceLung = getResizeImage(sliceLung)
img = np.hstack((sliceUint8.reshape((256, 256)), sliceLung.reshape((256, 256)), sliceLungPred.reshape((256, 256)))).astype("uint8")
#cv2.imshow("img", img)
#cv2.waitKey(0)
plt.imshow(img)
plt.gray()
plt.show()
# おまけ！
# Google ColaboratoryではGoogle Driveとは異なる独自のストレージが割り当てられ、
# 一定時間(12時間)経過すると、強制的にセッションが切れて、ダウンロードしたデータおよび学習モデル等が全て削除される(ノートの記述内容は残る)。
# なお、何もせず90分放置してもセッションが切れる。
# 学習したモデル等を保存したい場合、Google Driveに接続して、データを書き出すことも可能であるので、その例を示す。
# google driveマウントに必要なパッケージをインストール
!apt-get install -y -qq software-properties-common module-init-tools
!add-apt-repository -y ppa:alessandro-strada/ppa 2>&1 > /dev/null
!apt-get update -qq 2>&1 > /dev/null
!apt-get -y install -qq google-drive-ocamlfuse fuse
# Google colab用の認証トークン発行
from google.colab import auth
auth.authenticate_user()
#google-drive-ocamlfuseの許可取得
from oauth2client.client import GoogleCredentials
creds = GoogleCredentials.get_application_default()
import getpass
prompt = !google-drive-ocamlfuse -headless -id={creds.client_id} -secret={creds.client_secret} < /dev/null 2>&1 | grep URL
vcode = getpass.getpass(prompt[0] + '\n\nEnter verification code: ')
!echo {vcode} | google-drive-ocamlfuse -headless -id={creds.client_id} -secret={creds.client_secret}
# Googleドライブをマウントする
#マウントポイントを作成
!mkdir -p drive
!google-drive-ocamlfuse drive
print('Files in Drive:')
!ls drive/
#これ以降は残したいファイルはdriveへの保存やコピーが可能
#例1:(modelをdrive側のフォルダ内に保存するpythonコード)
# model.save_weights("drive/seminor/201902/model.h5")
#例2:(ファイルをdrive側のフォルダ内にコピーするコマンド)
#!cp semi_test.txt drive/seminor/201902/