In [ ]:
from __future__ import print_function
import os
from keras.layers import merge
from keras.layers.advanced_activations import LeakyReLU
from keras.optimizers import Nadam
from keras.metrics import MSE
from keras.layers.convolutional import Convolution2D, MaxPooling2D, ZeroPadding2D, AveragePooling2D
from keras.layers.core import Dense, Activation, Flatten
from keras.layers.normalization import BatchNormalization
from keras.models import Model
from keras.callbacks import ProgbarLogger, ModelCheckpoint
from keras.layers import Input
from keras.preprocessing.image import load_img, img_to_array, ImageDataGenerator
from sklearn.metrics import roc_auc_score, confusion_matrix, classification_report
from scipy.stats import spearmanr
import keras.backend as K
import numpy as np
import tensorflow as tf
import datetime
import scandir
import cv2
from PIL import Image
from ml_metrics import quadratic_weighted_kappa
os.environ['KERAS_BACKEND'] = 'tensorflow'
os.environ['CUDA_HOME'] = '/usr/local/cuda-7.5'
In [ ]:
def read_img(img_path):
'''This function returns a preprocessed image
'''
dim_ordering = K.image_dim_ordering()
img = load_img(img_path, target_size=(512, 512))
img = img_to_array(img, dim_ordering=dim_ordering)
if dim_ordering == 'th':
img = img[::-1, :, :]
else:
img = img[:, :, ::-1]
return img
In [ ]:
def identity_block(input_tensor, kernel_size, filters, stage, block, idx):
'''The identity_block is the block that has no conv layer at shortcut
# Arguments
input_tensor: input tensor
kernel_size: defualt 3, the kernel size of middle conv layer at main path
filters: list of integers, the nb_filters of 3 conv layer at main path
stage: integer, current stage label, used for generating layer names
block: 'a','b'..., current block label, used for generating layer names
'''
nb_filter1, nb_filter2, nb_filter3 = filters
if K.image_dim_ordering() == 'tf':
bn_axis = 3
else:
bn_axis = 1
conv_name_base = 'res' + str(stage) + block + '_branch' + '_' + 'idx'
bn_name_base = 'bn' + str(stage) + block + '_branch' + '_' + 'idx'
x = Convolution2D(nb_filter1, 1, 1, name=conv_name_base + '2a')(input_tensor)
x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2a')(x)
x = Activation('relu')(x)
x = Convolution2D(nb_filter2, kernel_size, kernel_size,
border_mode='same', name=conv_name_base + '2b')(x)
x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2b')(x)
x = Activation('relu')(x)
x = Convolution2D(nb_filter3, 1, 1, name=conv_name_base + '2c')(x)
x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2c')(x)
x = merge([x, input_tensor], mode='sum')
x = Activation('relu')(x)
return x
def conv_block(input_tensor, kernel_size, filters, stage, block, idx, strides=(2, 2)):
'''conv_block is the block that has a conv layer at shortcut
# Arguments
input_tensor: input tensor
kernel_size: defualt 3, the kernel size of middle conv layer at main path
filters: list of integers, the nb_filters of 3 conv layer at main path
stage: integer, current stage label, used for generating layer names
block: 'a','b'..., current block label, used for generating layer names
Note that from stage 3, the first conv layer at main path is with subsample=(2,2)
And the shortcut should have subsample=(2,2) as well
'''
nb_filter1, nb_filter2, nb_filter3 = filters
if K.image_dim_ordering() == 'tf':
bn_axis = 3
else:
bn_axis = 1
conv_name_base = 'res' + str(stage) + block + '_branch' + '_' + 'idx'
bn_name_base = 'bn' + str(stage) + block + '_branch' + '_' + 'idx'
x = Convolution2D(nb_filter1, 1, 1, subsample=strides,
name=conv_name_base + '2a')(input_tensor)
x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2a')(x)
x = Activation('relu')(x)
x = Convolution2D(nb_filter2, kernel_size, kernel_size, border_mode='same',
name=conv_name_base + '2b')(x)
x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2b')(x)
x = Activation('relu')(x)
x = Convolution2D(nb_filter3, 1, 1, name=conv_name_base + '2c')(x)
x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2c')(x)
shortcut = Convolution2D(nb_filter3, 1, 1, subsample=strides,
name=conv_name_base + '1')(input_tensor)
shortcut = BatchNormalization(axis=bn_axis, name=bn_name_base + '1')(shortcut)
x = merge([x, shortcut], mode='sum')
x = Activation('relu')(x)
return x
In [ ]:
# TODO: spatial parameters
def ResNet50(inp=None, idx=0):
'''Adapted from https://github.com/fchollet/deep-learning-models/blob/master/resnet50.py'''
# Determine proper input shape
if K.image_dim_ordering() == 'tf':
inp_shape=(512, 512, 3)
bn_axis = 3
else:
inp_shape=(3, 512, 512)
bn_axis = 1
dim_ordering = K.image_dim_ordering()
inp = Input(inp_shape)
with tf.device('/gpu:0'):
x = Convolution2D(64, 7, 7, subsample=(3, 3),
init='he_normal', border_mode='same', dim_ordering=dim_ordering,
name='conv1_{0}'.format(idx), input_shape=(inp_shape))(inp)
x = BatchNormalization(axis=bn_axis, name='bn_conv1_{0}'.format(idx))(x)
x = Activation('relu')(x) #LeakyReLU(alpha=0.5)
x = MaxPooling2D((3, 3), strides=(2, 2), dim_ordering=dim_ordering)(x)
x = conv_block(x, 3, [64, 64, 256], stage=2, block='a', idx=idx, strides=(1, 1))
x = identity_block(x, 3, [64, 64, 256], stage=2, block='b', idx=idx)
x = identity_block(x, 3, [64, 64, 256], stage=2, block='c', idx=idx)
with tf.device('/gpu:1'):
x = conv_block(x, 3, [128, 128, 512], stage=3, block='a', idx=idx)
x = identity_block(x, 3, [128, 128, 512], stage=3, block='b', idx=idx)
x = identity_block(x, 3, [128, 128, 512], stage=3, block='c', idx=idx)
x = identity_block(x, 3, [128, 128, 512], stage=3, block='d', idx=idx)
with tf.device('/gpu:2'):
x = conv_block(x, 3, [256, 256, 1024], stage=4, block='a', idx=idx)
x = identity_block(x, 3, [256, 256, 1024], stage=4, block='b', idx=idx)
x = identity_block(x, 3, [256, 256, 1024], stage=4, block='c', idx=idx)
x = identity_block(x, 3, [256, 256, 1024], stage=4, block='d', idx=idx)
x = identity_block(x, 3, [256, 256, 1024], stage=4, block='e', idx=idx)
x = identity_block(x, 3, [256, 256, 1024], stage=4, block='f', idx=idx)
with tf.device('/gpu:3'):
x = conv_block(x, 3, [512, 512, 2048], stage=5, block='a', idx=idx)
x = identity_block(x, 3, [512, 512, 2048], stage=5, block='b', idx=idx)
x = identity_block(x, 3, [512, 512, 2048], stage=5, block='c', idx=idx)
x = AveragePooling2D((7, 7), name='avg_pool_{0}'.format(idx))(x)
return Model(inp, x)
In [ ]:
patches_per_slice = 5
def get_mergenet():
if K.image_dim_ordering() == 'tf':
inp_shape=(512, 512, 3)
concat_axis = 3
else:
inp_shape=(3, 512, 512)
concat_axis = 1
patch_nets = []
input_list = []
for i in range(patches_per_slice):
inp = Input(inp_shape)
model = ResNet50(inp, i)
patch_nets.append(model)
input_list.append(inp)
out = merge(patch_nets, mode='ave')#, concat_axis=concat_axis)
out = Flatten()(out)
out = Dense(2, activation='softmax', name='fc')(out)
return Model(input_list, out)
In [ ]:
print('Started to build model at {0}'.format(datetime.datetime.utcnow()))
model = get_mergenet()
model.summary()
print('Finished building model at {0}'.format(datetime.datetime.utcnow()))
In [ ]:
#TODO: insert my own loss and metrics functions (one per output), equal to the ones specified by TUPAC
print('Started to compile model at {0}'.format(datetime.datetime.utcnow()))
model.compile(optimizer='adam',
loss='mse',
metrics=['mse'])
print('Finished compiling model at {0}'.format(datetime.datetime.utcnow()))
In [ ]:
def get_data(dir_path, res = 512):
annotations = open('../annotations/training_ground_truth.csv', 'r')
lines = annotations.readlines()
images = []
labels = []
for subdir, _, files in scandir.walk(dir_path):
subdir = subdir.replace('\\', '/') # windows path fix
subdir_split = subdir.split('/')
if len(subdir_split[3])>0: study_id = int(subdir_split[3].lstrip("0"))
else: continue
label = [float(lines[study_id].split(',')[0]), float(lines[study_id].split(',')[1])]
imgs = []
for f in files:
tiff_path = os.path.join(subdir, f)
if not tiff_path.endswith('.tiff'):
continue
img = read_img(tiff_path)
imgs.append(img)
images.append(imgs)
labels.append(label)
annotations.close()
images_list = []
for x in np.split(np.swapaxes(images, 0, 1), patches_per_slice):
images_list.append(np.squeeze(x))
return images_list, np.asarray(labels)
In [ ]:
# Load Data
images_train, labels_train = get_data('../example_images/train/')
images_val, labels_val = get_data('../example_images/val/')
print(images_train[0].shape, len(labels_train))
In [ ]:
# featurewise standardization, normalization and augmentation (horizontal and vertical flips)
def augment_data(X):
X = np.asarray(X).swapaxes(1, 0)
mean = np.mean(X, axis=0)
X -= mean
std = np.std(X, axis=0)
X /= (std + 1e-7)
X = X.swapaxes(0, 1)
if np.random.random() < 0.5:
X = flip_axis(X, 2)
if np.random.random() < 0.5:
X = flip_axis(X, 3)
return X
def flip_axis(X, axis):
X = np.asarray(X).swapaxes(axis, 0)
X = X[::-1, ...]
X = X.swapaxes(0, axis)
return X
print('Started data augmentation at {0}'.format(datetime.datetime.utcnow()))
images_train = augment_data(images_train)
images_val = augment_data(images_val)
print('Finished data augmentation at {0}'.format(datetime.datetime.utcnow()))
#print(np.mean(images_train[0]), np.std(images_train[0]))
In [ ]:
batch_size = 3
nb_epoch = 10
patches_per_slice = 5
#class_weights = {}
callbacks = [ProgbarLogger(),
ModelCheckpoint('mergenet_weights.hdf5', monitor='val_loss', verbose=1,
save_best_only=True, mode='auto')]
print('Started fitting at {0}'.format(datetime.datetime.utcnow()))
history = model.fit(images_train, labels_train, batch_size,
#samples_per_epoch=images_train.shape[0],
#class_weights = class_weights,
nb_epoch=nb_epoch, verbose=1,
validation_data=(images_val, labels_val),
#nb_val_samples = 1,
callbacks=callbacks)
print('Finished fitting at {0}'.format(datetime.datetime.utcnow()))
In [ ]:
# Decode predictions to one hot encoding
# Use ranking approach from https://github.com/ChenglongChen/Kaggle_CrowdFlower/blob/master/BlogPost/BlogPost.md
# The predicted values are mapped to classes based on the CDF of the ref values.
hist = np.bincount(labels_train[:,0].astype(int))
cdf = np.cumsum(hist) / float(sum(hist))
np.savetxt("../example_images/train/cdf.txt", cdf)
def getScore(pred, cdf, valid=False):
num = pred.shape[0]
output = np.asarray([4]*num, dtype=int)
rank = pred.argsort()
output[rank[:int(num*cdf[0]-1)]] = 1
output[rank[int(num*cdf[0]):int(num*cdf[1]-1)]] = 2
output[rank[int(num*cdf[1]):int(num*cdf[2]-1)]] = 3
if valid:
cutoff = [ pred[rank[int(num*cdf[i]-1)]] for i in range(3) ]
return output, cutoff
return output
In [ ]:
print('Metrics on validation set:\n------------------------------------------------------------\n')
pred = model.predict(images_val, verbose=0)
pred_class = getScore(pred[:,0], cdf)
print('Confusion Matrix:\n', confusion_matrix(labels_val[:,0], pred_class), '\n')
print(classification_report(labels_val[:,0], pred_class), '\n')
print("Quadratic Weighted Cohen's Kappa: ", quadratic_weighted_kappa(pred[:,0], labels_val[:,0]))
print("Spearman's Correlation Coëfficient: ", spearmanr(pred[:,1], labels_val[:,1])[0])
In [ ]:
print(history.history.keys())
In [ ]: