Import packages and functions¶
In [1]:
# Import packages
%matplotlib inline
from PIL import Image
import numpy as np
import os
import re
from skimage.color import gray2rgb
import matplotlib.pyplot as plt
from sklearn.utils import shuffle
!pip install tensorflow
!pip install keras
from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img
from keras.models import Sequential
from keras.layers import Dense, Dropout, Activation, Flatten, GaussianNoise, BatchNormalization, GlobalAveragePooling2D
from keras.layers import Conv2D, MaxPooling2D
from keras import Sequential
from keras.optimizers import Adam
import tensorflow as tf
from keras.backend.tensorflow_backend import set_session
from keras.preprocessing import image
from keras.models import Model
from keras import backend as K
from sklearn.metrics import confusion_matrix
from sklearn.metrics import roc_auc_score
!pip install git+https://github.com/raghakot/keras-vis.git --upgrade
from vis.visualization import visualize_cam, visualize_saliency, overlay
from keras import activations
from matplotlib import pyplot as plt
import matplotlib.cm as cm
import zipfile
from keras.models import model_from_json
import matplotlib as mpl
FIRST PART: DATA INGESTION¶
Import original images from local computer¶
These images come from a 4 year-old female with TSC. In total, there are 50 images consisting of 35 consecutive axial T2 MRI slices and 15 consecutive axial FLAIR MRI slices.
In [2]:
# Set the figure size
mpl.rcParams['figure.figsize'] = (16,10)
In [3]:
# Unzip files
with zipfile.ZipFile("TestCaseVT2.zip","r") as zip_ref:
zip_ref.extractall()
with zipfile.ZipFile("TestCaseVFLAIR.zip","r") as zip_ref:
zip_ref.extractall()
Path to original images folder¶
In [4]:
# Path to the folder with the original images
pathtoimagesT2test = './TestCaseVT2/'
pathtoimagesFLAIRtest = './TestCaseVFLAIR/'
SECOND PART: IMPORTATION OF FINAL DATA¶
In [5]:
# Functions to sort images with numbers within their name
def atoi(text):
return int(text) if text.isdigit() else text
def natural_keys(text):
return [ atoi(c) for c in re.split(r'(\d+)', text) ]
Import images and create labels for the T2 set¶
In [6]:
## T2
# Define the image size
image_size = (224, 224)
# Read in the test images for T2
T2test_images = []
T2test_dir = pathtoimagesT2test
T2test_files = os.listdir(T2test_dir)
T2test_files.sort(key=natural_keys)
# For each image
for f in T2test_files:
# Open the image
img = Image.open(T2test_dir + f)
# Resize the image so that it has a size 224x224
img = img.resize(image_size)
# Transform into a numpy array
img_arr = np.array(img)
# Transform from 224x224 to 224x224x3
if img_arr.shape == image_size:
img_arr = np.expand_dims(img_arr, 3)
img_arr = gray2rgb(img_arr[:, :, 0])
# Add the image to the array of images
T2test_images.append(img_arr)
# After having transformed all images, transform the list into a numpy array
T2test_X = np.array(T2test_images)
# Create an array of labels (as read by the radiologist)
T2test_y = np.array([[0], [0], [0], [0], [0], [1], [1], [1], [1], [1],
[1], [0], [0], [0], [0], [1], [1], [1], [1], [1],
[1], [0], [1], [1], [0], [0], [0], [0], [0], [0],
[0], [0], [0], [0], [0]])
# GPU expects values to be 32-bit floats
T2test_X = T2test_X.astype(np.float32)
# Rescale the values to be between 0 and 1
T2test_X /= 255.
In [7]:
T2test_X.shape
Out[7]:
In [8]:
# Example of an image to make sure they were converted right
plt.imshow(T2test_X[0])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
plt.show()
In [9]:
T2test_y.shape
Out[9]:
In [10]:
T2test_y[0]
Out[10]:
Import images and create labels for the FLAIR set¶
In [11]:
## FLAIR
# Define the image size
image_size = (224, 224)
# Read in the test images for FLAIR
FLAIRtest_images = []
FLAIRtest_dir = pathtoimagesFLAIRtest
FLAIRtest_files = os.listdir(FLAIRtest_dir)
FLAIRtest_files.sort(key=natural_keys)
# For each image
for f in FLAIRtest_files:
# Open the image
img = Image.open(FLAIRtest_dir + f)
# Resize the image so that it has a size 224x224
img = img.resize(image_size)
# Transform into a numpy array
img_arr = np.array(img)
# Transform from 224x224 to 224x224x3
if img_arr.shape == image_size:
img_arr = np.expand_dims(img_arr, 3)
img_arr = gray2rgb(img_arr[:, :, 0])
# Add the image to the array of images
FLAIRtest_images.append(img_arr)
# After having transformed all images, transform the list into a numpy array
FLAIRtest_X = np.array(FLAIRtest_images)
# Create an array of labels (as read by the radiologist)
FLAIRtest_y = np.array([[1], [1], [1], [1], [0], [0], [0], [0], [1], [1],
[0], [0], [0], [0], [0]])
# GPU expects values to be 32-bit floats
FLAIRtest_X = FLAIRtest_X.astype(np.float32)
# Rescale the values to be between 0 and 1
FLAIRtest_X /= 255.
In [12]:
FLAIRtest_X.shape
Out[12]:
In [13]:
# Example of an image to make sure they were converted right
plt.imshow(FLAIRtest_X[0])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
plt.show()
In [14]:
FLAIRtest_y.shape
Out[14]:
In [15]:
FLAIRtest_y[0]
Out[15]:
THIRD PART: VISUALIZE CLASS ACTIVATION MAPS AND SALIENCY MAPS¶
Load the model¶
In [16]:
# load model
json_file = open('InceptionV3.json', 'r')
loaded_model_json = json_file.read()
json_file.close()
model = model_from_json(loaded_model_json)
# load weights into new model
model.load_weights("InceptionV3.h5")
In [17]:
# Compile model
model.compile(optimizer = Adam(lr = 0.00025), loss = 'binary_crossentropy', metrics = ['accuracy'])
In [18]:
# Generate predictions on test data in the form of probabilities for T2
testInceptionV3T2 = model.predict(T2test_X, batch_size = 16)
testInceptionV3T2
Out[18]:
In [19]:
# Generate predictions on test data in the form of probabilities for FLAIR
testInceptionV3FLAIR = model.predict(FLAIRtest_X, batch_size = 16)
testInceptionV3FLAIR
Out[19]:
In [20]:
# Create the confusion matrix for T2
y_trueT2 = T2test_y
y_predInceptionV3T2 = testInceptionV3T2 > 0.5
confusion_matrix(y_trueT2, y_predInceptionV3T2, labels=[0,1])
Out[20]:
In [21]:
# Create the confusion matrix for FLAIR
y_trueFLAIR = FLAIRtest_y
y_predInceptionV3FLAIR = testInceptionV3FLAIR > 0.5
confusion_matrix(y_trueFLAIR, y_predInceptionV3FLAIR)
Out[21]:
In [22]:
# Calculate accuracy for T2
accuracy_InceptionV3T2 = (confusion_matrix(y_trueT2, y_predInceptionV3T2, labels=[0,1])[0, 0] + confusion_matrix(y_trueT2, y_predInceptionV3T2, labels=[0,1])[1, 1]) / (confusion_matrix(y_trueT2, y_predInceptionV3T2, labels=[0,1])[0, 0] + confusion_matrix(y_trueT2, y_predInceptionV3T2, labels=[0,1])[0, 1] + confusion_matrix(y_trueT2, y_predInceptionV3T2, labels=[0,1])[1, 0] + confusion_matrix(y_trueT2, y_predInceptionV3T2, labels=[0,1])[1, 1])
print('The accuracy in the test set is {}.'.format(accuracy_InceptionV3T2))
In [23]:
# Calculate accuracy for FLAIR
accuracy_InceptionV3FLAIR = (confusion_matrix(y_trueFLAIR, y_predInceptionV3FLAIR)[0, 0] + confusion_matrix(y_trueFLAIR, y_predInceptionV3FLAIR)[1, 1]) / (confusion_matrix(y_trueFLAIR, y_predInceptionV3FLAIR)[0, 0] + confusion_matrix(y_trueFLAIR, y_predInceptionV3FLAIR)[0, 1] + confusion_matrix(y_trueFLAIR, y_predInceptionV3FLAIR)[1, 0] + confusion_matrix(y_trueFLAIR, y_predInceptionV3FLAIR)[1, 1])
print('The accuracy in the test set is {}.'.format(accuracy_InceptionV3FLAIR))
Visualize the data¶
In [24]:
# Visualize the structure and layers of the model
model.layers
Out[24]:
In [25]:
# Iterate through the MRIs in T2
print('\n \n' + '\033[1m' + 'EACH ORIGINAL IMAGE IS ANALYZED WITH TWO METHODS: CLASS ACTIVATION MAP (UPPER ROW) AND SALIENCY MAP (LOWER ROW)' + '\033[0m' + '\n')
print('\033[1m' + 'FOR EACH METHOD, THE FIRST IMAGE IS THE ORIGINAL IMAGE, THE SECOND IMAGE IS THE MAP, AND THE THIRD IMAGE IS THE MAP SUPERIMPOSED ON THE ORIGINAL IMAGE WITH A TRANSPARENCY THAT IS PROPORTIONAL TO THE ESTIMATED PROBABILITY OF THE IMAGE HAVING TUBER(S) (HIGHER ESTIMATED PROBABILITIES PRODUCE CLEARLY SEEN MAPS OVERLAID ON THE ORIGINAL IMAGE AND LOWER ESTIMATED PROBABILITIES PRODUCE VERY TRANSPARENT MAPS OVERLAYED ON THE ORIGINAL IMAGE)' + '\033[0m'+ '\n \n \n \n')
for i in range(T2test_X.shape[0]):
# Print spaces to separate from the next image
print('\n \n \n \n \n \n \n \n')
# Print real classification of the image
print('\033[1m' + 'REAL CLASSIFICATION OF THE IMAGE: {}'.format('TUBER(S)' if y_trueT2[i][0]==1 else 'NO TUBER(S)') + '\033[0m')
# Print model classification and model probability of TSC
print('Model classification of this image: {} \nEstimated probability of tuber(s): {} \n'.format('TUBER(S)' if testInceptionV3T2[i][0]>0.5 else 'NO TUBER(S)', testInceptionV3T2[i][0]))
# Print title
print('\033[1m' + 'CLASS ACTIVATION MAP (UPPER ROW) AND SALIENCY MAP (LOWER ROW)' + '\033[0m')
# Original image
plt.subplot(2,3,1)
plt.imshow(T2test_X[i])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Heat map
plt.subplot(2,3,2)
heat_map = visualize_cam(model, layer_idx=300, filter_indices=None, seed_input=T2test_X[i])
plt.imshow(heat_map)
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Heat map superimposed on original image
plt.subplot(2,3,3)
plt.imshow(T2test_X[i])
plt.imshow(heat_map, alpha = 0.8 * testInceptionV3T2[i][0])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Original image
plt.subplot(2,3,4)
plt.imshow(T2test_X[i])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Heat map
heat_map = visualize_saliency(model, layer_idx=300, filter_indices=None, seed_input=T2test_X[i])
plt.subplot(2,3,5)
plt.imshow(heat_map)
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Heat map superimposed on original image
plt.subplot(2,3,6)
plt.imshow(T2test_X[i])
plt.imshow(heat_map, alpha = 0.8 * testInceptionV3T2[i][0])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Show the image and close it
plt.show()
plt.close()
SCROLL UP TO SEE THE GradCAM AND SALIENCY MAPS OF EACH IMAGE¶
This 4 year-old female had a medium tuber burden with tubers detected in 14 of 35 T2 MRI slices by our Neuroradiologist. The convolutional neural network classified images as having tuber(s) in 10 out of 35 T2 MRI slices and missed several T2 MRI slices with subtle tubers. Interestingly, in the missed tubers the convolutional neural network was focusing on the right are of tuber(s), but the estimated probability was too low to classify the MRI slice as a having tuber(s). Its accuracy in T2 was 0.57.¶
Each original image is analyzed with two methods: Gradient-weighted class activation maps (upper row) and saliency maps (lower row).¶
For each method, the first image is the original image, the second image is the map, and the third image is the map superimposed on the original image with a transparency that is proportional to the estimated probability of the image having tuber(s) (higher estimated probabilities produce clearly seen maps overlaid on the original image and lower estimated probabilities produce very transparent maps overlaid on the original image).¶
In [26]:
# Iterate through the MRIs in FLAIR
print('\n \n' + '\033[1m' + 'EACH ORIGINAL IMAGE IS ANALYZED WITH TWO METHODS: CLASS ACTIVATION MAP (UPPER ROW) AND SALIENCY MAP (LOWER ROW)' + '\033[0m' + '\n')
print('\033[1m' + 'FOR EACH METHOD, THE FIRST IMAGE IS THE ORIGINAL IMAGE, THE SECOND IMAGE IS THE MAP, AND THE THIRD IMAGE IS THE MAP SUPERIMPOSED ON THE ORIGINAL IMAGE WITH A TRANSPARENCY THAT IS PROPORTIONAL TO THE ESTIMATED PROBABILITY OF THE IMAGE HAVING TUBER(S) (HIGHER ESTIMATED PROBABILITIES PRODUCE CLEARLY SEEN MAPS OVERLAID ON THE ORIGINAL IMAGE AND LOWER ESTIMATED PROBABILITIES PRODUCE VERY TRANSPARENT MAPS OVERLAYED ON THE ORIGINAL IMAGE)' + '\033[0m'+ '\n \n \n \n')
for i in range(FLAIRtest_X.shape[0]):
# Print spaces to separate from the next image
print('\n \n \n \n \n \n \n \n')
# Print real classification of the image
print('\033[1m' + 'REAL CLASSIFICATION OF THE IMAGE: {}'.format('TUBER(S)' if y_trueFLAIR[i][0]==1 else 'NO TUBER(S)') + '\033[0m')
# Print model classification and model probability of TSC
print('Model classification of this image: {} \nEstimated probability of tuber(s): {} \n'.format('TUBER(S)' if testInceptionV3FLAIR[i][0]>0.5 else 'NO TUBER(S)', testInceptionV3FLAIR[i][0]))
# Print title
print('\033[1m' + 'CLASS ACTIVATION MAP (UPPER ROW) AND SALIENCY MAP (LOWER ROW)' + '\033[0m')
# Original image
plt.subplot(2,3,1)
plt.imshow(FLAIRtest_X[i])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Heat map
plt.subplot(2,3,2)
heat_map = visualize_cam(model, layer_idx=300, filter_indices=None, seed_input=FLAIRtest_X[i])
plt.imshow(heat_map)
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Heat map superimposed on original image
plt.subplot(2,3,3)
plt.imshow(FLAIRtest_X[i])
plt.imshow(heat_map, alpha = 0.8 * testInceptionV3FLAIR[i][0])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Original image
plt.subplot(2,3,4)
plt.imshow(FLAIRtest_X[i])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Heat map
heat_map = visualize_saliency(model, layer_idx=300, filter_indices=None, seed_input=FLAIRtest_X[i])
plt.subplot(2,3,5)
plt.imshow(heat_map)
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Heat map superimposed on original image
plt.subplot(2,3,6)
plt.imshow(FLAIRtest_X[i])
plt.imshow(heat_map, alpha = 0.8 * testInceptionV3FLAIR[i][0])
plt.grid(b=None)
plt.xticks([])
plt.yticks([])
# Show the image and close it
plt.show()
plt.close()
