In [ ]:
# coding=utf-8
from __future__ import print_function
import os
from keras.callbacks import ModelCheckpoint, EarlyStopping
from keras.optimizers import SGD
from sklearn.metrics import confusion_matrix
from scipy.stats import spearmanr
import openslide as ops
import cv2
import numpy as np
import datetime
import math
import json
from ml_metrics import quadratic_weighted_kappa
from keras.preprocessing.image import img_to_array
from keras.preprocessing.image import load_img
import keras.backend as K
import scandir
from sklearn import cross_validation
import matplotlib.pyplot as plt
from keras.layers import merge
from keras.layers.advanced_activations import LeakyReLU
from keras.layers.convolutional import Convolution2D
from keras.layers.convolutional import MaxPooling2D
from keras.layers.core import Dense, Activation, Flatten, MaxoutDense
from keras.layers.normalization import BatchNormalization
from keras.models import Model
from keras.layers import Input
import tensorflow as tf
os.environ['KERAS_BACKEND'] = 'tensorflow'
os.environ['CUDA_HOME'] = '/usr/local/cuda-7.5'
def read_img(img_path):
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
def splitAndSqueeze(images, patches_per_slice=3):
result = []
for img in np.split(images, patches_per_slice):
result.append(np.squeeze(img))
return result
# featurewise standardization, normalization and augmentation (horizontal and vertical flips)
def augment_data(X, patches_per_slice=3):
X = np.asarray(X).swapaxes(1, 0)
mean = np.mean(X, axis=0)
X = np.subtract(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 splitAndSqueeze(X)
def flip_axis(X, axis):
X = np.asarray(X).swapaxes(axis, 0)
X = X[::-1, ...]
X = X.swapaxes(0, axis)
return X
# 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.
def getScore(pred, cdf, valid=False):
num = pred.shape[0]
output = np.asarray([3]*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
def plotLearningCurve(history):
plt.plot(history.history['mean_squared_error'])
#plt.plot(history.history['val_mean_squared_error'])
plt.title('Learning Curve')
plt.ylabel('MSE')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.savefig('learning_curve_{0}.jpg'.format(datetime.datetime.utcnow()))
# Since we are looking for very local features, we should need a short amount of layers(2?), followed by a FCN
def get_patchblock(inp, idx):
if K.image_dim_ordering() == 'tf':
inp_shape=(512, 512, 3)
bn_axis = 3
else:
inp_shape=shape=(3, 512, 512)
bn_axis = 1
dim_ordering = K.image_dim_ordering()
out = Convolution2D(32, 7, 7, subsample=(2, 2),
init='he_normal', border_mode='same', dim_ordering=dim_ordering,
name='conv1_{0}'.format(idx), input_shape=(inp_shape))(inp)
out = BatchNormalization(axis=bn_axis, name='bn_conv1_{0}'.format(idx))(out)
out = Activation('relu')(out) #LeakyReLU(alpha=0.5)
out = MaxPooling2D((3, 3), strides=(2, 2), dim_ordering=dim_ordering)(out)
out = Convolution2D(32, 3, 3, subsample=(2, 2),
init='he_normal', border_mode='same', dim_ordering=dim_ordering,
name='conv2_{0}'.format(idx), input_shape=(inp_shape))(out)
out = BatchNormalization(axis=bn_axis, name='bn_conv2_{0}'.format(idx))(out)
out = Activation('relu')(out) #LeakyReLU(alpha=0.5)
out = Convolution2D(32, 3, 3, subsample=(2, 2),
init='he_normal', border_mode='same', dim_ordering=dim_ordering,
name='conv3_{0}'.format(idx), input_shape=(inp_shape))(out)
out = BatchNormalization(axis=bn_axis, name='bn_conv3_{0}'.format(idx))(out)
out = Activation('relu')(out) #LeakyReLU(alpha=0.5)
out = MaxPooling2D((3, 3), strides=(2, 2), dim_ordering=dim_ordering)(out)
out = Flatten()(out)
out = MaxoutDense(1, init='he_normal')(out)
return out
def get_mergenet(patches_per_slice=3):
print('Started to build model at {0}'.format(datetime.datetime.utcnow()))
if K.image_dim_ordering() == 'tf':
inp_shape=(512, 512, 3)
concat_axis = 3
else:
inp_shape=shape=(3, 512, 512)
concat_axis = 1
patch_nets = []
input_list = []
for i in range(patches_per_slice):
inp = Input(inp_shape)
model = get_patchblock(inp, i)
patch_nets.append(model)
input_list.append(inp)
out = merge(patch_nets, mode='ave')#, concat_axis=concat_axis)
out = MaxoutDense(1, init='he_normal')(out)
print('Finished building model at {0}'.format(datetime.datetime.utcnow()))
return Model(input_list, out)
In [ ]:
def get_data(dir_path, res = 512, limit=False):
annotations = open('../annotations/training_ground_truth.csv', 'r')
lines = annotations.readlines()
images_train = []
labels_train = []
images_val = []
labels_val = []
rs = cross_validation.ShuffleSplit(len(lines), n_iter = 1, test_size = 0.2, random_state = 17)
val_idx = []
count = 0
for train_index, val_index in rs:
val_idx.append(val_index)
for subdir, _, files in scandir.walk(dir_path):
for file in files:
study_id = int(file[9:12].lstrip("0"))
label = [float(lines[study_id-1].split(',')[0]), float(lines[study_id-1].split(',')[1])]
imgs = get_images(file[9:12])
if len(imgs)<3: continue
if (study_id in val_idx[0]):
images_val.append(np.asarray(imgs))
labels_val.append(np.asarray(label))
else:
images_train.append(np.asarray(imgs))
labels_train.append(np.asarray(label))
count += 1
if limit and count >=10:
break
if limit and count >=10:
break
annotations.close()
images_train = splitAndSqueeze(np.swapaxes(np.asarray(images_train), 0, 1))
images_val = splitAndSqueeze(np.swapaxes(np.asarray(images_val), 0, 1))
return images_train, images_val, np.asarray(labels_train), np.asarray(labels_val)
In [ ]:
def get_images(study_id, svs_dir='../../images/train/', roi_dir='../ROIs/'):
rois = open(roi_dir+'TUPAC-TR-{0}-ROI.csv'.format(study_id), 'r').readlines()
try:
slide = ops.OpenSlide(svs_dir+'{0}/{1}.svs'.format(study_id, study_id))
except ops.OpenSlideUnsupportedFormatError:
print('Not able to open svs file at: {0}'.format(svs_dir+'{0}/{1}.svs'.format(study_id, study_id)))
return []
imgs = []
for roi in rois:
roi = roi.split(',')
centroid = [int(roi[0])+int(math.ceil(int(roi[2])/2)), int(roi[1])+int(math.ceil(int(roi[3])/2))]
try:
level = 2
img = np.asarray(slide.read_region((centroid[0]-512, centroid[1]-512), 0, [1024, 1024]))
except ops.OpenSlideError:
print('Not able to open slide', study_id)
return []
res = 512
img = cv2.resize(img, (res, res), interpolation = cv2.INTER_AREA)
img = img[:, :, 0:3]
dim_ordering = K.image_dim_ordering()
img = img_to_array(img, dim_ordering=dim_ordering)
if dim_ordering == 'th':
img = img[::-1, :, :]
else:
img = img[:, :, ::-1]
imgs.append(img)
return imgs
In [ ]:
patches_per_slice = 3
model = get_mergenet()
model.summary()
In [ ]:
def spearmanloss(y_true, y_pred, batch_size=50):
length = batch_size
original_loss = -1*length*tf.reduce_sum(tf.mul(y_true, y_pred))-(tf.reduce_sum(y_true)*tf.reduce_sum(y_pred))
divisor = tf.sqrt((length * tf.reduce_sum(tf.square(y_true)) - tf.square(tf.reduce_sum(y_true))) *
(length * tf.reduce_sum(tf.square(y_pred)) - tf.square(tf.reduce_sum(y_pred)))
)
return tf.truediv(original_loss, divisor)
In [ ]:
opt = 'adam'#SGD(lr=0.1)
model.compile(optimizer=opt,loss='mse',metrics=['mse'])
images_train, images_val, labels_train, labels_val = get_data('../ROIs/')#, limit=True)
print(np.asarray(images_train).shape)
images_train = augment_data(images_train)
images_val = augment_data(images_val)
In [ ]:
batch_size = 50
nb_epoch = 30
callbacks = [ModelCheckpoint('mergenet_weights{0}.hdf5'.format(datetime.datetime.utcnow()), 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[:,1], batch_size,
nb_epoch=nb_epoch, verbose=1,
validation_data=(images_val, labels_val[:,1]),
callbacks=callbacks, shuffle=True)
print('Finished fitting at {0}'.format(datetime.datetime.utcnow()))
with open('01_mergenet_training_history{0}.json'.format(datetime.datetime.utcnow()), 'w') as f_out:
json.dump(history.history, f_out)
print('Metrics on validation set:\n------------------------------------------------------------\n')
pred = model.predict(images_val, verbose=0)
print("Spearman's Correlation Coëfficient: ", spearmanr(pred, labels_val[:,1])[0])
plotLearningCurve(history)
In [ ]:
%matplotlib inline
plt.hist(pred)
In [ ]:
plt.hist(labels_val[:,1])
In [ ]: