Publicado por In this project we develop a Deep Learning detector of Covid-19 in radiographs. For this purpose, we use images from the «Covid-chestxray-dataset» [3], generated by researchers from the Mila research group and the University of Montreal [4]. We also use images of radiographs of healthy and bacterial pneumonia patients extracted from Kaggle’s «Chest X-Ray Images (Pneumonia)» competition [5].
In total, we have a number of 426 images, divided into training (339 images), validation (42 images) and test (45 images) sets.
The partitions are divided into training (339 images), validation (42 images) and test (45 images).
The partitions are given in «.txt» lists, in which each image is assigned a tag:
- 0) Healthy
- 1) Covid-19
- 2) Pneumonia
Note: The results obtained by the models trained in this database are purely for educational purposes and cannot be used for actual diagnosis without clinical validation.
References¶
- María Climent, 2020 Covid-19: La Inteligencia Artificial De La Española Quibim Puede Acelerar El Diagnóstico Del Coronavirus
- Angel Alberich-bayarri,2020 Imagin, AI and Radiomix to understand and fight Coronavirus Covid-19
- Ieee8023/covid-chestxray-dataset
- Cohen, J.P., Morrison, P. and Dao, L., 2020. COVID-19 image data collection.
- Paul Mooney, 2019 Chest X-ray Images (pneumonia)
This notebook is run on Google Collab so we setup google drive to upload the images.
In [8]:
from google.colab import drive
drive.mount('/content/gdrive')
In [0]:
#Import libraries
import numpy as np
import re,shutil,os,timeit,glob
import matplotlib.pyplot as plt
import random
from IPython.display import Image
from sklearn.dummy import DummyClassifier
from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img
from keras.models import Sequential
from keras.layers import Dense, Conv2D, Flatten, Activation, Dropout, MaxPooling2D, BatchNormalization
from sklearn.metrics import accuracy_score
from keras.callbacks import ModelCheckpoint
from keras.optimizers import Adam,RMSprop,SGD,Adadelta
To load the images we will use the ImageDataGenerator from Keras, which generates images from the train, validation and test sets with the characteristics indicated.
To increase the train data we will generate moved or rotated images from the original dataset, also we will generate images with some noise. The images will be scaled to the appropriate size for each network architecture and the pixel values will be normalized.
In [0]:
#Function to clean txt files path.
def remChars(string):
#Remove first character
ret=string[2:]
ret=ret[:-1]
return ret
#Function to add noise to the images.
def add_noise(img):
VARIABILITY = 50
deviation = VARIABILITY*random.random()
noise = np.random.normal(0, deviation, img.shape)
img += noise
np.clip(img, 0., 255.)
return img
In [0]:
#Images path base
basepath="/content/gdrive/My Drive/"
In [0]:
#Let's create folder structure, train, test and validation with a folder per class inside them.
#Paths and txt files with the image names.
paths=['test','train','validation']
files=['testing.txt','training.txt','validation.txt']
#Read image names from the txt files
#Copy file in the proper folder.
for p,f in zip(paths, files):
file = open(basepath+f,"r")
imgfiles= file.readlines()
#Clean path and create folder structure
os.makedirs(basepath+p+"/COVID",exist_ok =True)
os.makedirs(basepath+p+"/HEALTHY",exist_ok =True)
os.makedirs(basepath+p+"/PNEUMONIA",exist_ok =True)
for s in imgfiles:
s=remChars(s)
if "COVID" in s:
shutil.copy(basepath+s, basepath+p+"/COVID")
elif "HEALTHY" in s:
shutil.copy(basepath+s,basepath+p+"/HEALTHY")
else:
shutil.copy(basepath+s,basepath+p+"/PNEUMONIA")
In [6]:
#Import images. Reduce size to 224x224
#Training dataset augmentation.
train_data_gen = ImageDataGenerator(
rescale=1./255,
rotation_range=25,
width_shift_range=0.3,
height_shift_range=0.3,
zoom_range=0.2,
horizontal_flip=True,
fill_mode='nearest',
preprocessing_function=add_noise)
#Only scale validation and test images.
validation_data_gen = ImageDataGenerator(rescale=1./255)
test_data_gen = ImageDataGenerator(rescale=1./255)
train_generator =train_data_gen.flow_from_directory(basepath+'train',
target_size=(224,224),
batch_size=32,
class_mode = 'categorical',
shuffle=True)
validation_generator = validation_data_gen.flow_from_directory(basepath+'validation',
target_size=(224,224),
batch_size=32,
class_mode = 'categorical',
shuffle=True)
test_generator = test_data_gen.flow_from_directory(basepath+'test',
target_size=(224,224),
batch_size=32,
class_mode = 'categorical',
shuffle=False)
In [0]:
#Show some transformed images examples.
img = load_img(train_generator.filepaths[0]) # this is a PIL image
x = img_to_array(img) # this is a Numpy array with shape (3, 150, 150)
x = x.reshape((1,) + x.shape) # this is a Numpy array with shape (1, 3, 150, 150)
i = 0
for batch in train_data_gen.flow(x, batch_size=1,
save_to_dir=basepath+'augmented', save_prefix='hi', save_format='jpeg'):
i += 1
if i > 2:
break
In [8]:
#Show images
for filename in glob.glob(basepath+'augmented/*.jpeg'): #assuming gif
display(Image(filename,width=150,height=150))
<
As the best option to classify images we will use convolutional networks.
For the classification of these images we will train models from different known CNN architectures (VGG,Inception,Densenet) and measure their accuracy and time performance to determine the most suitable one. We will use the pre-trained keras models, since our dataset is small and not similar to the dataset with which the network is trained, we will freeze some layers and re-train others.We will measure the accuracy of the results with the data from the testing dataset. We will also show the accuracy and loss plots to study the architectures.
We will first build a network with only one convolutional layer that we will use as a basis for defining further improvements.
The optimizers vary, as they have been adjusted in different test runs to try to get the most out of each type of network, the best results of each one are shown here.
In [0]:
#Simple model as baseline to measure performance improvements.
num_classes=3
model = Sequential()
model.add(Conv2D(11, kernel_size=3,strides=1, activation='relu', input_shape=(224,224,3)))
model.add(Flatten())
model.add((Dense(num_classes, activation='softmax')))
In [0]:
#Adam optimizer
optimizer=Adam(lr=2E-4)
# Compile the model
model.compile(optimizer=optimizer, loss='categorical_crossentropy', metrics=['accuracy'])
In [0]:
#Reset batch generator
train_generator.reset()
validation_generator.reset()
test_generator.reset()
In [14]:
#Run just 10 epochs and storing the execution time.
start_time = timeit.default_timer()
mfit=model.fit_generator(
train_generator,
epochs=15,
validation_data=validation_generator)
elapsedbaseline = timeit.default_timer() - start_time
In [15]:
#The elapsed time of this model is:
print('Baseline model elapsed time: '+str(elapsedbaseline) + ' seg.')
In [16]:
#Accuracy and Loss graphs for training and validation
# TODO
# summarize history for accuracy
plt.plot(mfit.history['accuracy'])
plt.plot(mfit.history['val_accuracy'])
plt.title('model accuracy')
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.show()
# TODO
# summarize history for loss
plt.plot(mfit.history['loss'])
plt.plot(mfit.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.show()
We see that the network is learning, albeit somewhat unevenly.
Let's check a VGG16 architecture, we import it from Keras.
In [20]:
from keras.applications.vgg16 import VGG16
vgg_c = VGG16(weights='imagenet', include_top=False,classes=3,input_shape=(224,224,3))
In [21]:
#Let's do transfer learning using the weights available in imagenet,
#having a small dataset, different from the original one with which the network was trained, we will run
#one part (only blocks 4 and 5).
for l in vgg_c.layers[:-8]:
l.trainable=False
for l in vgg_c.layers:
print(l.name+" "+str(l.trainable))
In [0]:
#Creamos las capas de salida
vgg_model=Sequential()
vgg_model.add(vgg_c)
vgg_model.add(Flatten())
vgg_model.add(Dense(4096,activation="relu"))
vgg_model.add(Dense(4096,activation="relu"))
vgg_model.add(Dense(3,activation="softmax"))
In [0]:
#Reset batches generators
train_generator.reset()
validation_generator.reset()
test_generator.reset()
In [24]:
#Let's define callbacks in order to save checkpoints just in case the process stops.
filepath = basepath+"vgg_model.h5"
checkpoint = ModelCheckpoint(filepath, monitor='loss', verbose=1, save_best_only=True, mode='min')
callbacks_list = [checkpoint]
#Optimizer setup
optimizer=Adam(lr=1E-6)
#Use same optimizer than in the base model
vgg_model.compile(loss='categorical_crossentropy', optimizer=optimizer, metrics=['accuracy'])
#Let's run 15 epochs, saving the elapsed time
start_time = timeit.default_timer()
mfitVGG=vgg_model.fit(
train_generator,
epochs=15,
validation_data=validation_generator,callbacks=callbacks_list)
elapsedVGG = timeit.default_timer() - start_time
In [25]:
#The elapsed time of this model has been:
print('VGG model elapsed time: '+str(elapsedVGG) + ' seg.')
In [27]:
#Train and validation accuarcy and loss graphs:
# TODO
# summarize history for accuracy
plt.plot(mfitVGG.history['accuracy'])
plt.plot(mfitVGG.history['val_accuracy'])
plt.title('model accuracy')
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.show()
# TODO
# summarize history for loss
plt.plot(mfitVGG.history['loss'])
plt.plot(mfitVGG.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.show()
With the VGG we see that we achieve greater accuracy in the validation dataset and that learning is performed more regularly.
We continue with an Inception architecture.
In [28]:
from keras.applications import InceptionV3
inception_c = InceptionV3(weights='imagenet', include_top=False,classes=3,input_shape=(299,299,3))
In [29]:
#We adapt the generators to the size of inception
train_generator =train_data_gen.flow_from_directory(basepath+'train',
target_size=(299,299),
batch_size=32,
class_mode = 'categorical',
shuffle=True)
validation_generator = validation_data_gen.flow_from_directory(basepath+'validation',
target_size=(299,299),
batch_size=32,
class_mode = 'categorical',
shuffle=True)
test_generator = test_data_gen.flow_from_directory(basepath+'test',
target_size=(299,299),
batch_size=32,
class_mode = 'categorical',
shuffle=False)
In [30]:
#Let's do transfer learning using the weights available in imagenet,
#having a small dataset, different from the original one with which the network was trained, we will run
#one part (only the last blocks).
for l in inception_c.layers[:-82]:
l.trainable=False
for l in inception_c.layers:
print(l.name+" "+str(l.trainable))
In [0]:
#We create the output layers
inception_model=Sequential()
inception_model.add(inception_c)
inception_model.add(Flatten())
inception_model.add(Dense(2048,activation="relu"))
inception_model.add(Dense(3,activation="softmax"))
In [0]:
#Reset batch generators
train_generator.reset()
validation_generator.reset()
test_generator.reset()
In [33]:
#We define callback to save checkpoints in case the process stops.
filepath = basepath+"inception.h5"
checkpoint = ModelCheckpoint(filepath, monitor='loss', verbose=1, save_best_only=True, mode='min')
callbacks_list = [checkpoint]
#Setup the optimizer
optimizer=Adam(lr=1E-6)
inception_model.compile(loss='categorical_crossentropy', optimizer=optimizer, metrics=['accuracy'])
#We run the model, saving the execution time, with 15 epochs
start_time = timeit.default_timer()
mfitxception=inception_model.fit(
train_generator,
epochs=15,
validation_data=validation_generator,callbacks=callbacks_list)
elapsedinception = timeit.default_timer() - start_time
In [34]:
elapsedinception
Out[34]:
In [36]:
#Accuracy and loss graphs for training and validation.
# TODO
# summarize history for accuracy
plt.plot(mfitxception.history['accuracy'])
plt.plot(mfitxception.history['val_accuracy'])
plt.title('model accuracy')
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.show()
# TODO
# summarize history for loss
plt.plot(mfitxception.history['loss'])
plt.plot(mfitxception.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.show()
The execution time in this case has increased considerably since it is the network with more layers. The accuracy is also gradually increasing, although the accuracy with the validation data does not reach the level of the VGG network.
We finally test with a DenseNet type network, in this case being a faster network, we run it completely, without using previously learned weights.
In [0]:
from keras.applications import DenseNet121
densenet = DenseNet121(weights=None, include_top=True,classes=3)
In [38]:
#for validation and test we just scale the image to the right size for the network
train_generator =train_data_gen.flow_from_directory(basepath+'train',
target_size=(224,224),
batch_size=32,
class_mode = 'categorical',
shuffle=True)
validation_generator = validation_data_gen.flow_from_directory(basepath+'validation',
target_size=(224,224),
batch_size=32,
class_mode = 'categorical',
shuffle=True)
test_generator = test_data_gen.flow_from_directory(basepath+'test',
target_size=(224,224),
batch_size=32,
class_mode = 'categorical',
shuffle=False)
In [0]:
train_generator.reset()
validation_generator.reset()
test_generator.reset()
In [46]:
#We define callback to save checkpoints in case the process stops.
filepath = basepath+"densenet.h5"
checkpoint = ModelCheckpoint(filepath, monitor='loss', verbose=1, save_best_only=True, mode='min')
callbacks_list = [checkpoint]
#Setup the optimizer
optimizer=SGD(lr=0.001)
densenet.compile(loss='categorical_crossentropy', optimizer=optimizer, metrics=['accuracy'])
#We run the model and save the runtime, with only 10 epochs
start_time = timeit.default_timer()
mfitdensenet=densenet.fit(
train_generator,
epochs=15,
validation_data=validation_generator,callbacks=callbacks_list)
densenetelapsed = timeit.default_timer() - start_time
In [47]:
densenetelapsed
Out[47]:
In [48]:
#Accuracy and loss graphs for train and validation data
# TODO
# summarize history for accuracy
plt.plot(mfitdensenet.history['accuracy'])
plt.plot(mfitdensenet.history['val_accuracy'])
plt.title('model accuracy')
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.show()
# TODO
# summarize history for loss
plt.plot(mfitdensenet.history['loss'])
plt.plot(mfitdensenet.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'validation'], loc='upper left')
plt.show()
In this case we see how we obtain a high accuracy in the validation dataset, but the learning does not seem very gradual, in fact in the training dataset there is hardly any improvement but it starts with a very high accuracy, while the accuracy in the validation does improve but with a lot of fluctuation.
Prediction and comparison¶
En esta sección se debe implementar la fase de test de los mejores modelos desarrollados anteriormente.
Let's make predictions on the test dataset to see the accuracy obtained with each of the networks with new data.
We start with the single-layer base model.
In [19]:
test_generator = test_data_gen.flow_from_directory(basepath+'test',
target_size=(224,224),
batch_size=32,
class_mode = 'categorical',
shuffle=False)
test_generator.reset()
#Prediction
prediction=model.predict_generator(test_generator)
#Cast to classes with probabilities
pclasses=[]
for p in prediction:
pclasses.append(np.argmax(p))
#We measure accuracy by comparing prediction with test data.
accuracy_score(test_generator.classes, pclasses)
Out[19]:
In this case the accuracy does not improve from 53%. We saw in the previous graphs that both the accuracy and loss of the network had remained somewhat stagnant even started to rise after the 15 stages, it is not expected that the network can improve much if more stages are added.
We continue with the VGG.
In [26]:
test_generator.reset()
# Realizamos la predicción con este modelo
prediction_vgg=vgg_model.predict(test_generator)
#Pasamos a clases con las probabilidades
pclasses_vgg=[]
for p in prediction_vgg:
pclasses_vgg.append(np.argmax(p))
#Medimos la precisión comparando la predicción con los datos de test
accuracy_score(test_generator.classes, pclasses_vgg)
Out[26]:
This model reaches a good accuracy, 84%, although at the cost of a considerable increase in processing time. We can see how we have been able to take good advantage of the previously learned weights provided by keras, from its training with imagenet.
Moreover, this network does seem to have prospects for improvement if we add more layers as both accuracy and loss improve regularly and gradually and do not seem to be at their maximum.
Let's go with Inception.
In [35]:
test_generator = test_data_gen.flow_from_directory(basepath+'test',
target_size=(299,299),
batch_size=32,
class_mode = 'categorical',
shuffle=False)
test_generator.reset()
# Let's go to the prediction with this model
prediction_inception=inception_model.predict(test_generator)
#Let's cast to classes with probabilities
pclasses_inception=[]
for p in prediction_inception:
pclasses_inception.append(np.argmax(p))
#We measure accuracy by comparing prediction with test data.
accuracy_score(test_generator.classes, pclasses_inception)
Out[35]:
Inception remains at 71% accuracy with the test data. It is not a bad figure and from the graphs it seems that it can continue to improve, although at the cost of a much longer processing time, given the amount of layers we are retraining.
Finally we performed predictions with DenseNet.
In [50]:
# Let's perform the prediction with this model
prediction_densenet=densenet.predict(test_generator)
#Let's cast to classes with probabilities
pclasses_densenet=[]
for p in prediction_densenet:
pclasses_densenet.append(np.argmax(p))
#We measure accuracy by comparing prediction with test data.
accuracy_score(test_generator.classes, pclasses_densenet)
Out[50]:
The result is quite good, 86%, better than VGG and with a lower processing time, although the graphs are a bit stranger since the network does not seem to learn in a stable way, since the training graphs are practically horizontal.
Discussion:
After evaluating the results of running the 4 networks, my conclusion is that for accuracy, run time and learning stability, the most suitable network for this dataset is VGG. We have been able to take good advantage of previously learned weights with other datasets and it seems that it can improve with more epochs. I think that by tuning the network more by modifying the number of retrained layers and adding epochs it can get more performance out of it than DenseNet. VGG is a deep network that uses small kernels to try to find more complex features in the network.
Explainability ¶
Let's try to visualize the important parts of one of the images for the CNN by implementing CAM (Class Activation Maps) with the library visualize_cam
Let's use a Densenet architecture.
In [34]:
#Gradcam algorithm implementation.
from keras.preprocessing.image import load_img, img_to_array
#Load and show the image
_img = load_img(basepath+"/train/COVID/ryct.2020200034.fig2.jpeg",target_size=(224,224))
plt.imshow(_img)
plt.show()
In [0]:
#We have to change scipy libraries in order to make visualize_cam to work.
!pip uninstall -y scipy
!pip install scipy==1.2.0
In [0]:
#Let's load the pretrained Densenet model.
from keras.applications.densenet import DenseNet121,preprocess_input
from keras.models import load_model
d_model=load_model(basepath+"densenet.h5")
In [50]:
#Let's generate a prediction with the image.
img = img_to_array(_img)
img = preprocess_input(img)
y_pred = d_model.predict(img[np.newaxis,...])
y_pred
class_idxs_sorted = np.argsort(y_pred.flatten())[::-1]
class_idxs_sorted
Out[50]:
In [54]:
#We have to set the linear activation in the last layer.
from vis.utils import utils
from keras.activations import linear
# Buscamos la última capa del modelo
layer_idx = utils.find_layer_idx(d_model, 'fc1000')
# Cambiamos su activación de softmax a linear
d_model.layers[layer_idx].activation = linear
d_model = utils.apply_modifications(d_model)
layer_idx
Out[54]:
In [58]:
#Let's search the last convolutional layer and its output.
penultima_capa = utils.find_layer_idx(d_model, "bn")
penultima_capa
Out[58]:
In [0]:
#No we can generate the heat map.
from vis.visualization import visualize_cam
class_idx = class_idxs_sorted[0]
seed_input = img
grad_top1 = visualize_cam(d_model, layer_idx, class_idx, seed_input,
penultimate_layer_idx = penultima_capa,#None,
backprop_modifier = None,
grad_modifier = None)
In [67]:
#Show the results.
def plot_map(grads):
fig, axes = plt.subplots(1,2,figsize=(14,5))
axes[0].imshow(_img)
axes[1].imshow(_img)
i = axes[1].imshow(grads,cmap="jet",alpha=0.8)
fig.colorbar(i)
plt.suptitle("heatmap")
plot_map(grad_top1)
