Deep Learning per la Diagnosi del Covid-19

Mentre scrivo questo post sono coricato a letto, con 38.7 di febbre, un mal di testa lancinante e le tonsille gonfie come due biglie. Adesso dovrei essere preoccupato, se non fosse per il fatto che le uniche interazioni sociali di presenza negli ultimi sei mesi le ho avute con il mio gatto Elon, quindi sto tranquillo e dedico questo tempo libero a quello che è il mio hobby preferito: binge watching su Github. Saltando da una repo all’altra ho trovato questo, una raccolta di radiografie del torace di pazienti aventi diverse malattie respiratorie, tra le quali anche il COVID-19. Forte dei miei mesi di full immersion nella Computer Vision per la realizzazione del mio ultimo video corso, ho deciso, come esercizio e come passatempo, di provare a realizzare una Rete Neurale per riconoscere il COVID-19 nelle radiografie.

premessa

Le Dipendenze

In questo progetto addestreremo una Rete Neurale Convoluzionale (ConvNet) e per farlo ci serviremo dei seguenti moduli Python:

• Tensorflow 2.0: popolarissima libreria per il Deep Learning realizzata da Google, la utilizzeremo per definire e addestrare la Rete Neurale.
• Scikit-learn: la più utilizzata libreria per il Machine Learning, la utilizzeremo per preprocessare le immagini.
• Pandas: libreria per l’analisi dati, la utilizzeremo per caricare e processare il file cvs con i metadati delle immagini.
• OpenCV: la più popolare libreria per la Computer Vision, la utilizzeremo per caricare e processare le immagini.
• Numpy: popolarissima libreria per il calcolo scientifico, la utilizzeremo per creare e manipolare array.
• Matplotlib e Seaborn: librerie per creare grafici e visualizzazioni con poche righe di codice, le utilizzeremo per visualizzare le radiografie e creare grafici con le prestazioni della rete.
• OS e Shutil: librerie per l’esecuzione di comandi di sistema, le utilizzeremo per leggere le immagini dalle directory e copiare quelle che ci interessano in nuove directory in maniera algoritmica.
• Math: modulo della standard library di Python che implementa diverse funzioni matematiche.

import os
import numpy as np
import pandas as pd

import cv2

from tensorflow.keras.preprocessing.image import ImageDataGenerator
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Conv2D, MaxPooling2D
from tensorflow.keras.layers import Dropout, Flatten, Dense

from sklearn.utils import shuffle
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import confusion_matrix

import matplotlib.pyplot as plt
import seaborn as sn

from shutil import copyfile

from math import sqrt

Definiamo un paio di costanti che ci serviranno durante il progetto.

BATCH_SIZE = 16 # Dimensione di ogni batch di immagini per l'algoritmo di ottimizzazione
IMG_SIZE = (200, 200) # Dimensione delle immagini

Il Dataset

Per addestrare la nostra Rete Neurale abbiamo bisogno di raccogliere un numero sufficiente di immagini, sia del caso positivo (pazienti affetti da COVID-19), che del caso negativo (pazienti sani). Creiamo la directory che conterrà le immagini.

if(not os.path.isdir("dataset")): # verifichiamo che la directory non esista già
  os.mkdir("dataset")
  os.mkdir("dataset/covid")
  os.mkdir("dataset/normal")

Radiografie di pazienti con il COVID-19

Iniziamo dalle immagini dei casi positivi. Cloniamo la repo con le immagini delle radiografie.

git clone https://github.com/ieee8023/covid-chestxray-dataset.git

All’interno della repo è presente un file chiamato metadata.csv, questo file contiene, diverse informazioni riguardo le immagini, carichiamolo in un Dataframe utilizzando Pandas e stampiamo le prime 5 righe.

df = pd.read_csv("covid-chestxray-dataset/metadata.csv") # creiamo il dataframe
df.head() # stampiamo le prime 5 righe

Dovremmo ottenere una tabella del genere:

Le informazioni che ci interessano sono:

• filename: contiene il nome del file .jpg con l’immagine
• finding: contiene il tipo di patologia diagnosticata
• view: contiene il tipo di proiezione della radiografia.

Filtriamo soltanto le righe dove la diagnosi è COVID e dove la proiezione è postero-anteriore (PA).

df = df[(df["finding"]=="COVID-19") & (df["view"]=="PA")]

Estraiamo all’interno di una lista i nomi dei vari file .jpg presenti nel dataframe filtrato, salviamo anche il numero di immagini in una variabile, ci servirà più avanti.

imgs_covid = list(df["filename"])
imgs_covid_count = len(imgs_covid)
imgs_covid_count

Nel momento in cui sto scrivendo questo articolo ci sono 99 immagini che possiamo utilizzare per il caso positivo.

Copiamo le varie immagini ottenute all’interno della cartella dataset/covid, che abbiamo creato in precedenza. Per farlo definiamo una funzione, che riutilizzeremo più avanti anche per le immagini del caso negativo.

def copy_samples(imgs_list, imgs_path, cls):
  for img in imgs_list:
    copyfile(imgs_path+img, "dataset/"+cls+"/"+img)

copy_samples(imgs_covid, "covid-chestxray-dataset/images/", "covid")

RADIOGRAFIE DI PAZIENTI SANI

Passiamo alle radiografie di pazienti sani, per ottenerle utilizzeremo questo dataset presente su Kaggle.

Possiamo scaricare il dataset direttamente dalla sua pagina su Kaggle, oppure utilizzando le API. Io scelgo questa seconda opzione, in modo da rendere il notebook facilmente eseguibile per tutti su Google Colaboratory. Se non lo abbiamo già, utilizziamo pip per installare il modulo kaggle. Quindi da terminale:

pip install kaggle

Creiamo il file kaggle.json con le credenziali per l’utilizzo delle API, per ottenere le tue credenziali guarda qui.

user = "IL_TUO_NOME_UTENTE" 
key = "LA_TUA_API_KEY"

if '.kaggle' not in os.listdir('/root'):
    !mkdir ~/.kaggle
!touch /root/.kaggle/kaggle.json
!chmod 666 /root/.kaggle/kaggle.json
with open('/root/.kaggle/kaggle.json', 'w') as f:
    f.write('{"username":"%s","key":"%s"}' % (user, key))
!chmod 600 /root/.kaggle/kaggle.json

Ora usiamo il modulo kaggle per scaricare il dataset, sempre da terminale:

kaggle datasets download -d paultimothymooney/chest-xray-pneumonia

Tramite il programma unzip (o una qualsiasi alternativa per Windows) estraiamo l’archivio che abbiamo scaricato.

unzip -q chest-xray-pneumonia.zip

All’interno dell’archivio sono presenti tre cartelle, train, val e test. Dentro queste cartelle le immagini sono a loro volta contenute in altre due cartelle, NORMAL e PNEUMONIA. Le immagini che a noi interessano sono quelle di pazienti sani, quindi quelle contenute all’interno della cartella NORMAL. Per creare un dataset bilanciato dobbiamo selezionare un numero di immagini pari quello delle immagini del caso positivo (questo è il motivo per il quale abbiamo salvato tale numero in precedenza). Quindi:

1. Leggiamo i nomi di tutte le immagini nella cartella train/NORMAL/
2. Mescoliamo la lista così ottenuta
3. Selezioniamo dalla lista mescolata solo i primi N elementi, dove N è il numero di immagini che abbiamo a disposizione per il caso positivo.

from random import shuffle

imgs_normal = os.listdir("chest_xray/chest_xray/train/NORMAL")
shuffle(imgs_normal)
imgs_normal = imgs_normal[:imgs_covid_count]

Ora utilizziamo la funzione definita in precedenza per copiare le immagini nella directory del nostro dataset.

copy_samples(imgs_normal, "chest_xray/train/NORMAL/", "normal")

Preprocessing delle Immagini

Abbiamo il dataset ! Prima di poter addestrare la rete neurale dobbiamo preparare le immagini. Iniziamo caricando immagini e labels in due array numpy in un loop, all’interno del quale eseguiremo anche le seguenti trasformazioni:

• conversione in grayscale (bianco e nero): trattandosi di radiografie, i colori non ci interessano.
• ridimensionamento delle immagini con la dimensione che abbiamo definito nella variable IMG_SIZE.
• normalizzazione delle immagini portando tutti i pixel in un range di valori che va da 0 a 1.
• equalizzazione dell’istogramma: in modo da limitare differenze di contrasto/luminosità nelle immagini.
• codifica dei label assegnando al caso normale il valore 0 (negativo) e al caso covid il valore 1 (positivo).

X = []
y = []

encoding = [("normal",0),("covid",1)] # definiamo la codifica

for folder, label in encoding:
  for img_name in os.listdir("dataset/"+folder):
    img = cv2.imread("dataset/"+folder+"/"+img_name, cv2.IMREAD_GRAYSCALE) # carichiamo l'immagine in un array numpy
    img = cv2.equalizeHist(img) # equalizziamo l'istogramma
    img = cv2.resize(img, IMG_SIZE)/255. # ridimensionamento e normalizzazione
    X.append(img)
    y.append(label)

X = np.array(X)
y = np.array(y)

Per esser sicuri di aver fatto tutto correttamente, visualizziamo su schermo qualche immagine. Definiamo una funzione per farlo.

def show_samples(X):

  fig = plt.figure() # creiamo una nuova figura

  for i in range(X.shape[0]):
    plot = fig.add_subplot(1,X.shape[0],i+1) # aggiugniamo un subplot per l'immagine
    plt.imshow(X[i]) # mostriamo l'immagine
    plt.axis("off") # rimuoviamo i valori sulle assi

Ora utilizziamo la funzione per mostrare 5 radiografie del caso negativo…

show_samples(X[1:6])

e 5 del caso positivo…

show_samples(X[-6:-1])

Gli esempi all’interno degli array sono ordinati, il primo 50% contiene esempi negativi e il restante 50% contiene esempi positivi, questo ci ha aiutato a stampare le immagini dei due casi, ma non va assolutamente bene per addestrare una rete neurale, in quanto questa potrebbe incappare in dei cicli. Ora potremmo mescolare gli esempi presenti nei due array sfruttando la funzione shuffle di sklearn, ma possiamo anche farne a meno, dato che dobbiamo dividere gli array per addestramento e test utilizzando la funzione train_test_split e questa mescola gli array in automatico.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, stratify=y)
print(X_train.shape[0])
print(X_test.shape[0])

Per completare, aggiungiamo un’altra dimensione agli array pari al numero di canali delle immagini, cioè 1 dato che si tratta di immagini in bianco e nero, questo va fatto perché Tensorflow richiede l’input della rete in questo formato.

X_train = X_train.reshape(X_train.shape[0], X_train.shape[1], X_train.shape[2], 1)
X_test = X_test.reshape(X_test.shape[0], X_test.shape[1], X_test.shape[2], 1)

DATA AUGMENTATION

Il numero di immagini raccolte nella repository su Github aumenta di giorno in giorno, però allo stato attuale sono ancora troppo poche per poter addestrare un qualsiasi modello di Machine Learning, figuriamoci una Deep Convolutional Neural Network. Per provare ad aggirare questa limitazione utilizzeremo due semplici tecniche di data augmentation:

• generiamo nuove immagini ruotando quelle esistenti di un massimo di 15 gradi, in senso orario e antiorario.
• generiamo nuove immagini aumentando/riducendo la luminosità di quelle esistenti.

Per farlo sfruttiamo i generatori di Keras.

datagen = ImageDataGenerator(
        rotation_range=15,
        fill_mode="nearest"
)

train_generator = datagen.flow(
        X_train,
        y_train,
        batch_size=BATCH_SIZE,)

Adesso possiamo addestrare il modello sul generatore di immagini.

Addestramento della Rete Neurale

Finalmente siamo arrivati alla parte più eccitante del progetto, l’addestramento della Rete Neurale Convoluzionale.

Per prima cosa definiamone l’architettura, dopo qualche tentativo questa è quella che mi ha portato ad un risultato ottimale:

• Un primo strato convoluzionale con 32 filtri
• Uno strato di Pooling
• Un secondo strato convoluzionale con 32 filtri
• Uno strato di Pooling
• Uno strato denso con 64 nodi
• Infine lo strato di output con un unico nodo.

Avendo davvero pochi dati il rischio di overfitting è molto alto.

L’overfitting è una condizione in cui un modello di machine learning memorizza la struttura dei dati, piuttosto che utilizzarla per apprendere, e quindi fallisce nel generalizzare su dati non visti durante la fase di addestramento.

Per contrastare l’overfitting utilizziamo una tecnica chiamata Dropout, che consiste nel disattivare in maniera casuale un percentuale dei nodi di uno strato ad ogni iterazione dell’algoritmo di ottimizzazione, limitando situazioni di co-adaptations (perdonami l’inglesismo ma non ho idea di quale sia il termine italiano), cioè situazioni in cui nodi tendono a farsi carico degli errori di altri nodi portando ad una condizione di overfitting.

model = Sequential()
model.add(Conv2D(filters=32, kernel_size=3, activation="relu", input_shape=(IMG_SIZE[0], IMG_SIZE[1], 3)))
model.add(MaxPooling2D(pool_size=2))
model.add(Dropout(0.6))

model.add(Conv2D(filters=32, kernel_size=3, activation="relu"))
model.add(MaxPooling2D(pool_size=2))
model.add(Dropout(0.6))

model.add(Flatten())
model.add(Dense(64, activation="relu"))
model.add(Dropout(0.6))
model.add(Dense(1, activation="sigmoid"))

Siamo quasi pronti, per ultimare la rete neurale dobbiamo configurare il processo di apprendimento:

• trattandosi di un problema di classificazione binaria, come funzione di costo utilizziamo la binary crossentropy (o log loss).
• come algoritmo di ottimizzazione utilizziamo l’Adam, una variante del Gradient Descent che utilizza un learning rate adattivo e sfrutta il Momentum.
• aggiungiamo l’accuracy tra le metriche da visualizzare durante l’addestramento

model.compile(loss='binary_crossentropy',
              optimizer='adam',
              metrics=['accuracy'])

Adesso è tutto pronto per poter avviare l’addestramento, facciamolo per 100 epoche.

model.fit(X_train, y_train, epochs=100)

Questo è il mio output per le ultime 10 epoche:

L’addestramento si è concluso e le metriche sul generatore, log loss e accuracy,  sono ottime, ma per essere sicuri della qualità del nostro modello dobbiamo verificarle sul set di test.

VALUTIAMO IL MODELLO

Valutiamo il modello calcolando accuracy log loss sia sul set di addestramento che sul set di test, possiamo farlo automaticamente utilizzando il metodo evaluate.

metrics_train = model.evaluate(X_train, y_train)
metrics_test = model.evaluate(X_test, y_test)

print("Train Accuracy = %.4f - Train Loss = %.4f" % (metrics_train[1], metrics_train[0]))
print("Test Accuracy = %.4f - Test Loss = %.4f" % (metrics_test[1], metrics_test[0]))

Queste sono le mie metriche:

Il risultato è ottimo (anche troppo forse) ma c’è un ma ! Quando si tratta di diagnosi mediche, il classificare un caso positivo come negativo è molto più grave che classificare un caso negativo come positivo, in quanto in quest’ultimo caso il paziente potrebbe essere messo in quarantena senza un reale motivo, fino a quando ulteriori esami non accerteranno che è sano,  mentre nel primo caso al paziente verrebbero negate le giuste cure e verrebbe rimandato a casa, dove, in questo caso, potrebbe anche infettare altre persone. Quale ti sembra lo scenario peggiore ?

Per vedere dove la nostra rete neurale sta sbagliando creiamo e visualizziamo una matrice di confusione.

y_pred = model.predict_classes(X_test)
cm = confusion_matrix(y_pred, y_test)
df_cm = pd.DataFrame(cm, index = ["Predicted Normal","Predicted Covid"],
                  columns = ["Normal","Covid"])
sn.heatmap(df_cm, annot=True)

Questa è l’immagine che viene generata:

La matrice di confusione è così formata:

• Primo riquadro (in alto a sinistra): casi negativi classificati come negativi, cioè i pazienti sani classificati come sani.
• Secondo riquadro (in alto a destra): casi positivi classificati come negativi, cioè i pazienti affetti dal Covid-19 classificati come sani.
• Terzo riquadrato (in basso a sinistra): casi negativi classificati come positivi, cioè i pazienti sani classificati come affetti dal Covid-19.
• Quarto riquadrato (in basso a destra):  casi positivi classificati come positivi, cioè pazienti con il Covid-19 classificati come affetti dal Covid-19.

Il secondo riquadro è quello alla quale dobbiamo prestare massima attenzione, come vedi è presente un falso negativo, il che è un male per le motivazioni che ho spiegato sopra. Una semplice tecnica che possiamo adottare è aumentare la sensibilità della rete, ad esempio classificando come positive le immagini che forniscono una probabilità maggiore del 45%, piuttosto che del 50% che si è soliti adottare. Proviamo.

y_proba = model.predict(X_test)
y_pred = np.where(y_proba>0.45,1,0)
cm = confusion_matrix(y_pred, y_test)

df_cm = pd.DataFrame(cm, index = ["Predicted Normal","Predicted Covid"],
                  columns = ["Normal","Covid"])
sn.heatmap(df_cm, annot=True)

Ed ecco l’immagine:

Testiamo la Rete Neurale

Andiamo a testare la nostra rete neurale sul campo, eseguendo delle predizioni sulle radiografie. Purtroppo non abbiamo altre immagini a disposizione quindi eseguiremo le predizioni sempre sul set di test.

pred_proba = model.predict(X_test)
print(np.round(pred_proba,3))

Ecco qui l’output:

In questo caso abbiamo ottenuto le predizioni sotto forma di probabilità di appartenenza alla classe positiva.

Facciamo una cosa migliore, stampiamo le immagini delle radiografie all’interno di una griglia, scrivendo nella parte alta di ogni immagine la diagnosi della nostra rete neurale.

labels = ["normal", "covid"]
colors = ["green","red"]

pred = model.predict_classes(X_test)

fig = plt.figure(figsize=(14,16))
fig.subplots_adjust(wspace=0, hspace=0)
n_cols = int(sqrt(X_test.shape[0]))
n_rows = X_test.shape[0]//n_cols+1

for i in range(X_test.shape[0]):

  plot = fig.add_subplot(n_rows, n_cols, i+1)
  plt.imshow(X_test[i])

  class_index = pred[i][0]
  plt.text(5, 25, labels[class_index], fontsize=14, color=colors[class_index], fontdict={'weight': 'bold'})
  plt.axis("off")

Questo è il risultato finale:

CONSIDERAZIONI FINALI

• I risultati che abbiamo ottenuto sono fantastici, anche troppo fantastici e nel machine learning non bisogna mai farsi coinvolgere da risultati emozionanti perché l’inghippo è sempre dietro l’angolo. In questo caso, avendo davvero poche immagini a disposizione, non possiamo esser sicuri di cosa ha realmente appreso la rete, anche perché la difficoltà di interpretabilità delle Reti Neurali è uno dei principali limit del Deep Learning. Comunque come ho già detto all’inizio dell’articolo il mio scopo non voleva essere quello di creare uno strumento diagnostico all’avanguardia per il COVID-19, onestamente non penso di essere in grado di farlo, bensì mostrarti come il Deep Learning potrebbe essere applicato in questi casi.
• Se l’articolo ti è piaciuto fammelo sapere con un commento, se sarete in molti ne realizzerò una continuazione in cui tenterò di debuggare la rete per capire cosa realmente sta apprendendo.
• Se ti interessa l’argomento Deep Learning/Computer Vision forse il mio ultimo videocorso potrebbe fare al caso tuo, da pochi giorni abbiamo avviato l’iniziativa IORESTOACASA, se acquisti un corso con il codice promozionale IORESTOACASA e lo completi entro il 30 Aprile il successivo te lo regalo io :).