在上一次肺炎X光片的預測中,我們通過神經網絡來識別患者胸部的X光片，用於檢測患者是否患有肺炎。這是一個典型的神經網絡圖像分類在醫學領域中的運用。
另外，神經網絡的圖像分割在醫學領域中也有着很重要的用作。接下來，我們要演示如何在氣胸患者的X光片上，分割出氣胸患者患病區的部位和形狀。
那么就讓我們來正式開始了。
第一步：導入需要的 Python 包
import sys
import cv2
import pydicom
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import patches as patches
from glob import glob
from tqdm import tqdm
from plotly.offline import download_plotlyjs, init_notebook_mode, iplot
from plotly import subplots
from plotly.graph_objs import *
from plotly.graph_objs.layout import Margin, YAxis, XAxis
init_notebook_mode()
import tensorflow as tf
from tensorflow import reduce_sum
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Input, Conv2D, MaxPool2D, UpSampling2D, Concatenate, Flatten, Add
from tensorflow.keras.losses import binary_crossentropy
from tensorflow.keras import callbacks
from sklearn.model_selection import train_test_split
# 數據增強庫
from albumentations import (
Compose, HorizontalFlip, CLAHE, HueSaturationValue,
RandomBrightness, RandomContrast, RandomGamma,OneOf,
ToFloat, ShiftScaleRotate,GridDistortion, ElasticTransform, JpegCompression, HueSaturationValue,
RGBShift, RandomBrightness, RandomContrast, Blur, MotionBlur, MedianBlur, GaussNoise,CenterCrop,
IAAAdditiveGaussianNoise,GaussNoise,OpticalDistortion,RandomSizedCrop
)
# 設置使用90%的顯存，避免顯存OOM錯誤
config = tf.compat.v1.ConfigProto()
config.gpu_options.per_process_gpu_memory_fraction = 0.9
session = tf.compat.v1.Session(config=config)
%matplotlib inline
第二步：導入數據庫
圖像分割的數據一共分為兩部分：
訓練用的圖片
圖片中需要分割的部分，稱之為 mask
這次我們訓練的數據是以 DCM 文件存儲。
DCM 是一種數位成像，廣泛運用於醫學領域，但並不局限於醫學，DCM 本身是一種特殊的圖像文件，它可以用來存儲各種圖像信息。
DCM 文件是遵循 DICOM 標准的一種文件
我們的 mask 部分的數據存儲在 csv 文件中，csv 文件大家都比較熟悉, 這里就不做介紹了。
2.1 導入 mask 數據
首先我們來看下存放 mask 數據的 csv 文件中的 mask 數據
使用 pandas 的 read_csv 接口讀取 train-rle.csv 文件。
我們先查看其中的頭部5條數據，其中我們可以看到這個 csv 文件中存放了兩列，一列是 ImageId ， 一列是 EncodedPixels 。
ImageId 這一列比較好理解，是訓練數據的 id，對應的是 dcm 文件的文件名。
rles_df = pd.read_csv('pneumothorax-segmentation/train-rle.csv')
rles_df.columns = ['ImageId', 'EncodedPixels']
rles_df.head()
這里要對 EncodedPixels 這一列做下說明：
EncodedPixels 實際存放的就是 mask 的像素數據，這些像素數據是以 RLE 編碼存放的。
接下來，我們需要定義一個函數來將 RLE 編碼的數據還原成 mask 圖片數據。
def rle2mask(rle, width, height):
mask= np.zeros(width* height)
array = np.asarray([int(x) for x in rle.split()])
starts = array[0::2]
lengths = array[1::2]
current_position = 0
for index, start in enumerate(starts):
current_position += start
mask[current_position:current_position+lengths[index]] = 255
current_position += lengths[index]
return mask.reshape(width, height)
2.2 導入 DCM 文件
接下來我們將DCM讀入並存儲到字典中，方便以后查看跟使用。我們還將之前讀入的 mask 數據也合並到相應的 ImageId 的字典中。
在訓練數據中，如果胸片沒有被 mask 標記，表示這個病例他並不患有氣胸。通過 EncodedPixels 中的數據，將是否是氣胸的患者記錄到 has_pneumothorax 這一字段中。
def dicom_to_dict(dicom_data, file_path, rles_df, encoded_pixels=True):
data = {}
# Parse fields with meaningful information
data['patient_name'] = dicom_data.PatientName
data['patient_id'] = dicom_data.PatientID
data['patient_age'] = int(dicom_data.PatientAge)
data['patient_sex'] = dicom_data.PatientSex
data['pixel_spacing'] = dicom_data.PixelSpacing
data['file_path'] = file_path
data['id'] = dicom_data.SOPInstanceUID
# look for annotation if enabled (train set)
if encoded_pixels:
encoded_pixels_list = rles_df[rles_df['ImageId']==dicom_data.SOPInstanceUID]['EncodedPixels'].values
pneumothorax = False
for encoded_pixels in encoded_pixels_list:
if encoded_pixels != ' -1':
pneumothorax = True
data['encoded_pixels_list'] = encoded_pixels_list
data['has_pneumothorax'] = pneumothorax
data['encoded_pixels_count'] = len(encoded_pixels_list)
return data
train_fns = sorted(glob('pneumothorax-segmentation/dicom-images-train/*/*/*.dcm'))
train_metadata_df = pd.DataFrame()
train_metadata_list = []
for file_path in tqdm(train_fns):
dicom_data = pydicom.dcmread(file_path)
train_metadata = dicom_to_dict(dicom_data, file_path, rles_df)
train_metadata_list.append(train_metadata)
train_metadata_df = pd.DataFrame(train_metadata_list)
train_metadata_df.head()
第三步：數據可視化
我們在讀取完數據以后，接下來就進行數據情況的查看。
3.1 隨機挑選病例樣本
我們隨機挑選了幾個病例。我們在每個病例上打出了年齡，性別以及是否是氣胸患者。
num_img = 4
subplot_count = 0
fig, ax = plt.subplots(nrows=1, ncols=num_img, sharey=True, figsize=(num_img*10,10))
for index, row in train_metadata_df.sample(n=num_img).iterrows():
dataset = pydicom.dcmread(row['file_path'])
ax[subplot_count].imshow(dataset.pixel_array, cmap=plt.cm.bone)
# label the x-ray with information about the patient
ax[subplot_count].text(0,0,'Age:{}, Sex: {}, Pneumothorax: {}'.format(row['patient_age'],row['patient_sex'],row['has_pneumothorax']),
size=26,color='white', backgroundcolor='black')
subplot_count += 1
結果如圖示：
我們在看下 mask 圖像在相對應的病例中的位置：
我們分三組來顯示
第一組我們將原始胸片圖像中用紅色的框框出 mask 的最小包圍盒. 然后將mask 數據部分用不同的顏色區分
第二組我們將原始圖像做直方圖均衡化處理，讓胸片對比度更加清晰。
第三組我們直接顯示原始圖像
通過觀察我們看到，如果沒有一定的專業知識，根本無法區分跟看出氣胸的具體位置。
def bounding_box(img):
# return max and min of a mask to draw bounding box
rows = np.any(img, axis=1)
cols = np.any(img, axis=0)
rmin, rmax = np.where(rows)[0][[0, -1]]
cmin, cmax = np.where(cols)[0][[0, -1]]
return rmin, rmax, cmin, cmax
def plot_with_mask_and_bbox(file_path, mask_encoded_list, figsize=(20,10)):
pixel_array = pydicom.dcmread(file_path).pixel_array
clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(16, 16))
clahe_pixel_array = clahe.apply(pixel_array)
# use the masking function to decode RLE
mask_decoded_list = [rle2mask(mask_encoded, 1024, 1024).T for mask_encoded in mask_encoded_list]
fig, ax = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(20,10))
# print out the xray
ax[0].imshow(pixel_array, cmap=plt.cm.bone)
# print the bounding box
for mask_decoded in mask_decoded_list:
# print out the annotated area
ax[0].imshow(mask_decoded, alpha=0.3, cmap="Reds")
rmin, rmax, cmin, cmax = bounding_box(mask_decoded)
bbox = patches.Rectangle((cmin,rmin),cmax-cmin,rmax-rmin,linewidth=1,edgecolor='r',facecolor='none')
ax[0].add_patch(bbox)
ax[0].set_title('With Mask')
結果如圖示：
3.2 氣胸患者的數據分布
接下來，我們需要查看患有氣胸的數據和未患有氣胸的數據的分布情況。
nok_count = train_metadata_df['has_pneumothorax'].sum()
ok_count = len(train_metadata_df) - nok_count
x = ['No Pneumothorax','Pneumothorax']
y = [ok_count, nok_count]
trace0 = Bar(x=x, y=y, name = 'Ok vs Not OK')
nok_encoded_pixels_count = train_metadata_df[train_metadata_df['has_pneumothorax']==1]['encoded_pixels_count'].values
trace1 = Histogram(x=nok_encoded_pixels_count, name='# of annotations')
fig = subplots.make_subplots(rows=1, cols=2)
fig.append_trace(trace0, 1, 1)
fig.append_trace(trace1, 1, 2)
fig['layout'].update(height=400, width=900, title='Pneumothorax Instances')
iplot(fig)
結果如圖示：
3.3 氣胸患者的年齡分布
在讓我們通過年齡的角度來看下氣胸患者的分布情況。
train_male_df = train_metadata_df[train_metadata_df['patient_sex']=='M']
train_female_df = train_metadata_df[train_metadata_df['patient_sex']=='F']
pneumo_pat_age = train_metadata_df[train_metadata_df['has_pneumothorax']==1]['patient_age'].values
no_pneumo_pat_age = train_metadata_df[train_metadata_df['has_pneumothorax']==0]['patient_age'].values
pneumothorax = Histogram(x=pneumo_pat_age, name='has pneumothorax')
no_pneumothorax = Histogram(x=no_pneumo_pat_age, name='no pneumothorax')
fig = subplots.make_subplots(rows=1, cols=2)
fig.append_trace(pneumothorax, 1, 1)
fig.append_trace(no_pneumothorax, 1, 2)
fig['layout'].update(height=400, width=900, title='Patient Age Histogram')
iplot(fig)
結果如圖示：
3.4 氣胸患者的性別分布
讓我們通過性別的角度來查看下氣胸患者的分布情況。
train_male_df = train_metadata_df[train_metadata_df['patient_sex']=='M']
train_female_df = train_metadata_df[train_metadata_df['patient_sex']=='F']
male_ok_count = len(train_male_df[train_male_df['has_pneumothorax']==0])
female_ok_count = len(train_female_df[train_female_df['has_pneumothorax']==0])
male_nok_count = len(train_male_df[train_male_df['has_pneumothorax']==1])
female_nok_count = len(train_female_df[train_female_df['has_pneumothorax']==1])
ok = Bar(x=['male', 'female'], y=[male_ok_count, female_ok_count], name='no pneumothorax')
nok = Bar(x=['male', 'female'], y=[male_nok_count, female_nok_count], name='has pneumothorax')
data = [ok, nok]
layout = Layout(barmode='stack', height=400)
fig = Figure(data=data, layout=layout)
iplot(fig, filename='stacked-bar')
結果如圖示：
m_pneumo_labels = ['no pneumothorax','has pneumothorax']
f_pneumo_labels = ['no pneumothorax','has pneumothorax']
m_pneumo_values = [male_ok_count, male_nok_count]
f_pneumo_values = [female_ok_count, female_nok_count]
colors = ['#FEBFB3', '#E1396C']
fig = {
"data": [
{
"values": m_pneumo_values,
"labels": m_pneumo_labels,
"domain": {"column": 0},
"name": "Male",
"hoverinfo":"label+percent",
"hole": .4,
"type": "pie"
},
{
"values": f_pneumo_values,
"labels": f_pneumo_labels,
"textposition":"inside",
"domain": {"column": 1},
"name": "Female",
"hoverinfo":"label+percent",
"hole": .4,
"type": "pie"
}],
結果如圖示：
數據可視化是分析數據的一種重要手段，通過上面幾個例子，給大家展示了一些比較常用的數據可視化的方法。
第四步：數據清洗
下面我們來看下，我們的數據內是否含有無效的數據，無效的數據指的是我們胸片圖片跟 mask 上不一致，也可以說我們的胸片並未被標記標簽。
在第二步驟中有說明 mask 是以 RLE 編碼的，如果是氣胸患者，那么他的 RLE 數據段是有值的，如果他不是我們是以 -1 來標示。
在剛才讀取 dcm 的函數里，我們把 EncodedPixels 字段的數據長度給記錄下來，正常的數據長度必須是 >0 。
因此我們可以簡單的查看下記錄 EncodedPixels 長度的 encoded_pixels_count 字段中值是否為 0，來簡單的過濾下我們的非正常數據。
missing_vals = train_metadata_df[train_metadata_df['encoded_pixels_count']==0]['encoded_pixels_count'].count()
print("Number of x-rays with missing masks: {}".format(missing_vals))
我們可以看到，有37份數據是沒有標簽，后面我們需要刪除它們。
第五步：准備訓練數據
我們先來定義一些我們將要使用到的一些常數。
# 圖像大小
img_size = 256
# batch size
batch_size = 16
# 卷積kernel的大小
k_size = 3
# 訓練數據跟驗證數據的分割比例
val_size = .1
5.1 數據生成類
class DataGenerator(tf.keras.utils.Sequence):
def __init__(self, file_path_list, labels, batch_size=32,
img_size=256, channels=1, shuffle=True, augmentations=None):
self.file_path_list = file_path_list
self.labels = labels
self.batch_size = batch_size
self.img_size = img_size
self.channels = channels
self.shuffle = shuffle
self.augment = augmentations
self.on_epoch_end()
def __len__(self):
return int(np.floor(len(self.file_path_list)) / self.batch_size)
def __getitem__(self, index):
indexes = self.indexes[index*self.batch_size:(index+1)*self.batch_size]
file_path_list_temp = [self.file_path_list[k] for k in indexes]
X, y = self.__data_generation(file_path_list_temp)
if self.augment is None:
return X,np.array(y)/255
else:
im,mask = [],[]
for x,y in zip(X,y):
augmented = self.augment(image=x, mask=y)
im.append(augmented['image'])
mask.append(augmented['mask'])
return np.array(im),np.array(mask)/255
5.2 數據增強
數據增強有助於提高數據的數量，因為你每做一次變換相當於得到了一張新的圖片；
同時也能提高模型的泛化能力，因為你的數據分布相較於沒做數據增強的數據的分布更加的廣泛。
AUGMENTATIONS_TRAIN = Compose([
HorizontalFlip(p=0.5),
OneOf([
RandomContrast(),
RandomGamma(),
RandomBrightness(),
], p=0.3),
OneOf([
ElasticTransform(alpha=120, sigma=120 * 0.05, alpha_affine=120 * 0.03),
GridDistortion(),
OpticalDistortion(distort_limit=2, shift_limit=0.5),
], p=0.3),
RandomSizedCrop(min_max_height=(156, 256), height=img_size, width=img_size, p=0.25),
ToFloat(max_value=1)
],p=1)
AUGMENTATIONS_VALIDATION = Compose([
ToFloat(max_value=1)
],p=1)
第六步：定義模型
在圖像分割中, 有很多模型都可以達到很好的效果. 今天我們選用的模型是 UNet 的變形。
使用tensorflow的keras接口定義我們的 ResUNet 。
'batch normalization layer with an optinal activation layer'
x = tf.keras.layers.BatchNormalization()(x)
if act == True:
x = tf.keras.layers.Activation('relu')(x)
return x
def conv_block(x, filters, kernel_size=3, padding='same', strides=1):
'convolutional layer which always uses the batch normalization layer'
conv = bn_act(x)
conv = Conv2D(filters, kernel_size, padding=padding, strides=strides)(conv)
return conv
def stem(x, filters, kernel_size=3, padding='same', strides=1):
conv = Conv2D(filters, kernel_size, padding=padding, strides=strides)(x)
conv = conv_block(conv, filters, kernel_size, padding, strides)
shortcut = Conv2D(filters, kernel_size=1, padding=padding, strides=strides)(x)
shortcut = bn_act(shortcut, act=False)
output = Add()([conv, shortcut])
return output
def residual_block(x, filters, kernel_size=3, padding='same', strides=1):
res = conv_block(x, filters, k_size, padding, strides)
res = conv_block(res, filters, k_size, padding, 1)
shortcut = Conv2D(filters, kernel_size, padding=padding, strides=strides)(x)
shortcut = bn_act(shortcut, act=False)
output = Add()([shortcut, res])
return output
定義我們的模型對象，並且 Summary 一下看下模型的細節。
model = ResUNet(img_size)
model.compile(optimizer="adam", loss=bce_dice_loss, metrics=[iou_metric])
model.summary()
第七步：開始訓練
epochs=70
callback = LearningRateCallbackBuilder(nb_epochs=epochs,nb_snapshots=1,init_lr=1e-3)
history = model.fit_generator(generator=training_generator, validation_data=validation_generator, callbacks=callback.get_callbacks(), epochs=epochs, verbose=2)
第八步：訓練結果查看
# 模型的IoU
plt.figure(figsize=(16,4))
plt.subplot(1,2,1)
plt.plot(history.history['iou_metric'][1:])
plt.plot(history.history['val_iou_metric'][1:])
# 模型的loss
plt.subplot(1,2,2)
plt.plot(history.history['loss'][1:])
plt.plot(history.history['val_loss'][1:])
plt.ylabel('val_loss')
訓練數據的 loss 跟驗證數據的 loss 的走勢對比，訓練數據的 IoU 跟驗證數據的 IoU 的走勢對比。
通過上面兩個走勢我們可以看出，我們的模型在不斷的收斂的。
count = 0
for i in range(0,30):
if count <= 15:
x, y = validation_generator.__getitem__(i)
predictions = model.predict(x)
for idx, val in enumerate(x):
if y[idx].sum() > 0 and count <= 15:
img = np.reshape(x[idx]* 255, (img_size, img_size))
mask = np.reshape(y[idx]* 255, (img_size, img_size))
pred = np.reshape(predictions[idx], (img_size, img_size))
pred = pred > 0.5
pred = pred * 255
plot_train(img, mask, pred)
count += 1
通過上述圖片，我們可以看到氣胸的陰影面積和位置，已經被分離出來了。但是，某些參數還需要進一步的調整。
大家可以登陸矩池雲國內領先的GPU雲共享平台，選擇demo鏡像，進行該氣胸分割案例嘗試。