U-Netを用いた胸部X線画像からの肺野セグメンテーション

##環境

* Ubuntu 18.04
* Anaconda 1.9.6
* Tensorflow-gpu 1.12.0
* Keras-gpu 2.2.4
* GeForce GTX 1050 Ti(totalMemory: 3.95GiB freeMemory: 3.89GiB) + AkiTiONODE

※GPUを使いましたが、CPUでもできます。(確かめてませんが。)

##必要なもの

日本放射線技術学会「標準ディジタル画像データベース[胸部腫瘤陰影像]」Segmentation01
http://imgcom.jsrt.or.jp/minijsrtdb/

![Screenshot from 2019-05-04 06-22-26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379410/58ac1c9f-e3ed-ad50-9fe1-a268e6544371.png)

画像サイズ　256*256、色深度　Gray 8 bitのデータ(一番上のもの)を選択してダウンロード。
ダウンロードしたら解凍。ここではわかりやすくデスクトップへ。
Segmentation01.zipをデスクトップに移動して、右クリックから「ここで展開」を実行。

![Screenshot from 2019-05-04 06-25-29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379410/e73c06d8-cc68-bc39-f3a0-b8d5fcdbce57.png)

Segmentation01フォルダの中身は、とても親切にトレーニングとテストデータが分かれており、かつ、それぞれのフォルダの中身は、org(original image), label(mask image)に整理されている。作成・提供者に感謝。

![Screenshot from 2019-05-04 06-26-42.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379410/f70fb96a-286d-a550-68c6-b51ff06a59e0.png)


##コード

少し長いです。。

```
import os
import sys
import random
import warnings

import numpy as np
import pandas as pd
import glob

import matplotlib.pyplot as plt

from skimage.io import imread, imshow

from keras.models import Model, load_model
from keras.layers import Input
from keras.layers.core import Lambda
from keras.layers import BatchNormalization
from keras.layers.convolutional import Conv2D, Conv2DTranspose
from keras.layers.pooling import MaxPooling2D
from keras.layers.merge import concatenate
from keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau
from keras import backend as K
from keras.optimizers import Nadam
from keras.layers.advanced_activations import LeakyReLU

import tensorflow as tf

# いくつかパラメータをセット
IMG_WIDTH = 256　#今回利用する画像のサイズ
IMG_HEIGHT = 256　#今回利用する画像のサイズ
IMG_CHANNELS = 1　#今回利用する画像のチャンネル（グレースケールなので１）
# ご自身の環境に合わせて変更して下さい。
# 例えば、Win10の方は、'C:\Users\あなたのユーザ名\Desktop\Segmentation01\train'などになります。
TRAIN_PATH = '/home/tatsunidas/デスクトップ/Segmentation01/train/'
TEST_PATH = '/home/tatsunidas/デスクトップ/Segmentation01/test/'
# UserWarningを無視する設定（なくてもOK）
warnings.filterwarnings('ignore', category=UserWarning, module='skimage')

# オリジナルデータとマスクデータのパスを取得、テストデータも。
train_imgs = glob.glob(TRAIN_PATH+'org/*.png')
train_masks = glob.glob(TRAIN_PATH+'label/*.png')
test_imgs = glob.glob(TEST_PATH+'org/*.png')
test_masks = glob.glob(TEST_PATH+'label/*.png')

# オリジナルとマスク画像の枚数が違っていたらコメントを出力
num_of_train_imgs = len(train_imgs)
num_of_train_masks = len(train_masks)
if num_of_train_imgs != num_of_train_masks:
    print('invalid datasets, please check train data')
# 各画像を配列で取得しておく
# image
X_train = np.zeros((num_of_train_imgs, IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS), dtype=np.uint8)
# mask
Y_train = np.zeros((num_of_train_masks, IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS), dtype=np.bool)
print('Getting and resizing train images and masks ... ')
sys.stdout.flush()
for n in range(num_of_train_imgs):
    X_train[n] = imread(train_imgs[n],as_gray=True).reshape(IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS)
    Y_train[n] = imread(train_masks[n],as_gray=True).reshape(IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS)

# テスト画像を配列で取得しておく
X_test = np.zeros((len(test_imgs), IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS), dtype=np.uint8)
Y_test = np.zeros((len(test_masks), IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS), dtype=np.uint8)
print('Getting test images ... ')
sys.stdout.flush()
for n in range(len(test_imgs)):
    X_test[n] = imread(test_imgs[n],as_gray=True).reshape(IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS)
    Y_test[n] = imread(test_masks[n],as_gray=True).reshape(IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS)
print('Preparing is Done!')

# 念のため画像のペアを確認（何回も実行する場合はコメントアウトしてください）
ix = random.randint(0, len(train_imgs))
imshow(X_train[ix].reshape(IMG_HEIGHT, IMG_WIDTH))
plt.show()
imshow(np.squeeze(Y_train[ix].reshape(IMG_HEIGHT, IMG_WIDTH)))
plt.show()

# mean IoUという評価指標を利用する
# Define IoU metric
def mean_iou(y_true, y_pred):
    prec = []
    for t in np.arange(0.5, 1.0, 0.05):
        y_pred_ = tf.to_int32(y_pred > t)
        score, up_opt = tf.metrics.mean_iou(y_true, y_pred_, 2)
        K.get_session().run(tf.local_variables_initializer())
        with tf.control_dependencies([up_opt]):
            score = tf.identity(score)
        prec.append(score)
    return K.mean(K.stack(prec), axis=0)

# U-Net modelをビルド
# 少し手を加えています。騙されないように注意して下さい笑（悪意はありません）
inputs = Input((IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS))
s = Lambda(lambda x: x / 255.) (inputs)

c1 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (s)
c1 = BatchNormalization() (c1)
c1 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c1)
p1 = MaxPooling2D((2, 2)) (c1)

c2 = Conv2D(32, (3, 3), kernel_initializer='he_normal', padding='same') (p1)
c2 = LeakyReLU()(c2)
c2 = BatchNormalization() (c2)
c2 = Conv2D(32, (3, 3), kernel_initializer='he_normal', padding='same') (c2)
c2 = LeakyReLU()(c2)
p2 = MaxPooling2D((2, 2)) (c2)

c3 = Conv2D(64, (3, 3), kernel_initializer='he_normal', padding='same') (p2)
c3 = LeakyReLU()(c3)
c3 = BatchNormalization() (c3)
c3 = Conv2D(64, (3, 3), kernel_initializer='he_normal', padding='same') (c3)
c3 = LeakyReLU()(c3)
p3 = MaxPooling2D((2, 2)) (c3)

c4 = Conv2D(128, (3, 3), kernel_initializer='he_normal', padding='same') (p3)
c4 = LeakyReLU()(c4)
c4 = BatchNormalization() (c4)
c4 = Conv2D(128, (3, 3), kernel_initializer='he_normal', padding='same') (c4)
c4 = LeakyReLU()(c4)
p4 = MaxPooling2D(pool_size=(2, 2)) (c4)

c5 = Conv2D(256, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (p4)
c5 = BatchNormalization() (c5)
c5 = Conv2D(256, (3, 3), activation='elu',kernel_initializer='he_normal', padding='same') (c5)

u6 = Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same') (c5)
u6 = concatenate([u6, c4])
c6 = Conv2D(128, (3, 3), kernel_initializer='he_normal', padding='same') (u6)
c6 = LeakyReLU()(c6)
c6 = BatchNormalization() (c6)
c6 = Conv2D(128, (3, 3), kernel_initializer='he_normal', padding='same') (c6)
c6 = LeakyReLU()(c6)

u7 = Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same') (c6)
u7 = concatenate([u7, c3])
c7 = Conv2D(64, (3, 3), kernel_initializer='he_normal', padding='same') (u7)
c7 = LeakyReLU()(c7)
c7 = BatchNormalization() (c7)
c7 = Conv2D(64, (3, 3), kernel_initializer='he_normal', padding='same') (c7)
c7 = LeakyReLU()(c7)

u8 = Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same') (c7)
u8 = concatenate([u8, c2])
c8 = Conv2D(32, (3, 3), kernel_initializer='he_normal', padding='same') (u8)
c8 = LeakyReLU()(c8)
c8 = BatchNormalization() (c8)
c8 = Conv2D(32, (3, 3), kernel_initializer='he_normal', padding='same') (c8)
c8 = LeakyReLU()(c8)

u9 = Conv2DTranspose(16, (2, 2), strides=(2, 2), padding='same') (c8)
u9 = concatenate([u9, c1], axis=3)
c9 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (u9)
c9 = BatchNormalization(momentum=0.9) (c9)
c9 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c9)

outputs = Conv2D(1, (1, 1), activation='sigmoid') (c9)

model = Model(inputs=[inputs], outputs=[outputs])
# loss関数は、Dice Coefficientを使う場合もあるが、今回はbinary_crossentropy
# metricsに'binary_accuracy',mean_iouを設定して学習精度を観察する
model.compile(optimizer=Nadam(lr=2e-5), loss='binary_crossentropy', metrics=['binary_accuracy',mean_iou])
model.summary()

# モデルのトレーニング
# 損失値が小さくならなくなったら学習をストップするかどうか
# earlystopper = EarlyStopping(patience=10, verbose=1)
# 一番いい状態のときにモデルを保存する設定
checkpointer = ModelCheckpoint('model-u-net.h5', verbose=1, save_best_only=True)
# 損失値が小さくならなくなったら学習率を小さくする設定
reduceLr = ReduceLROnPlateau(monitor='val_loss',factor=0.1,patience=7)
# 上記設定を加味してトレーニング
results = model.fit(X_train, 
                    Y_train, 
                    validation_split=0.1, #バリデーションデータをトレーニングデータから自動分割
                    batch_size=3, 
                    epochs=50, 
                    callbacks=[reduceLr, checkpointer])

# 学習後にテストデータで予測結果を確認
model = load_model('model-u-net.h5', custom_objects={'mean_iou': mean_iou})
preds_test = model.predict(X_test, verbose=1)
# ピクセルワイズの評価結果を0.5の閾値で0 or 255に変換（評価結果をバイナリ画像に変換）。
preds_test_t = (preds_test > 0.5).astype(np.uint8)

# 結果をランダムに１つ取得して表示
ix = random.randint(0, len(preds_test_t))
imshow(X_test[ix].reshape(IMG_HEIGHT, IMG_WIDTH))
plt.show()
imshow(Y_test[ix].reshape(IMG_HEIGHT, IMG_WIDTH))
plt.show()
imshow(np.squeeze(preds_test_t[ix].reshape(IMG_HEIGHT, IMG_WIDTH)))
plt.show()
```

##結果

1s 32ms/step - loss: 0.0539 - binary_accuracy: 0.9843 - mean_iou: 0.8347 - val_loss: 0.0837 - val_binary_accuracy: 0.9730 - val_mean_iou: 0.8358

オリジナル
![Figure_11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379410/b9e299c3-0741-30d9-63c7-04579e5d29e9.png)
教師マスク
![Figure_mask11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379410/1960c490-8629-e871-7a74-45a0cb56f09d.png)
今回の評価サンプル
![Figure_11result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379410/e6e69ad5-0a3e-897e-00b7-5811522d55f9.png)

もう少し良く出来そうな予感。他にもDeepResU-Netなども試してみたい。

##この記事を書いたモチベーション

備忘録として。勉強として。

##その他

コメント歓迎です。誤りなど修正致します。

##参考URL
コードの殆どの部分：https://www.kaggle.com/keegil/keras-u-net-starter-lb-0-277
Res U-Netを参考にしたい:https://github.com/DuFanXin/deep_residual_unet