Deep learning e segmentazione per la biologia cellulare

Nell’era del digitale anche in medicina ogni giorno sono prodotte numerose immagini e sebbene la loro classificazione risulta ampliamente studiata la segmentazione e il riconoscimento di oggetti sono meno approfonditi. Con il recente sviluppo dell’intelligenza artificiale e degli algoritmi delle reti neurali convoluzionali è possibile ottenere ottimi risultati. In questo articolo vengono analizzate le differenti reti neurali in Matlab e la possibilità di riconoscere cellule in microscopia ottica. Viene poi introdotta anche la possibilità di segmentare bordo e zona interna estendono l’algoritmo dalle cellule in microscopia al riconoscimento di nei per la dermatoscopia.

Nella biologia è presente da diversi anni la necessità di automatizzare il riconoscimento cellulare. La possibilità di riconoscere una cellula precisa all’interno di una qualsivoglia immagine in microscopia è un’abilità molto richiesta. Per decenni tale compito è stato svolto da operatori umani ma con l’aumento esponenziale della potenza di calcolo dei moderni computer è stato possibile introdurre diversi algoritmi. Oggi giorno esistono diverse tecniche per identificare oggetti all’interno delle immagini, rientrando nel campo della segmentazione, ma la maggior parte di esse risentono di problemi legati alla qualità, forma, dimensione, convessità e molte altre proprietà geometriche dell’oggetto di interesse. L’ultimo sviluppo vede l’applicazione delle reti neurali e del transfer learning per automatizzare completamente la segmentazione ottenendo prestazioni straordinarie.

Background

Segmentazione

La segmentazione è formalmente definita come “dividere un’immagine in un set di regioni non sovrapponibili che unite danno l’intera immagine”. Ovvero è un insieme di tecniche per veicolare ed estrarre l’informazione dalle immagini attraverso alcune proprietà, tra cui quelle morfologiche assegnado ad ogni pixel un’etichetta univoca, associata ad un insieme di classi. Esistono diversi metodi e principi di segmentazione e in tutti si parla di foreground e background. Il foreground è proprio ciò che risulta di interesse e quindi lo si vuole separare dal background. Esistono quattro categorie principali di metodi:

• Pixel-based, basati sulla luminanza dei singoli pixel
• Edge-based, sfruttano bordi e criteri di discontinuità
• Region-based, utilizzano criteri di somiglianza tra regioni vicine
• Model-based, richiedono un modello geometrico dell’oggetto da cercare

A questi si aggiungono i metodi supervisionati che richiedono reti neurali e transfer learning per automatizzare il riconoscimento.

Metriche di valutazione

Oltre all’immagine di partenza è necessario avere un ground truth (GT) ovvero una maschera di verità che esprime quale pixel dell’immagine è effettivamente un pixel dell’oggetto di interesse. Sulla base di questa maschera e dell’identificazione che segue la segmentazione vengono definite alcune metriche di qualità.

Ci saranno dei pixel falsi positivi (FP) ovvero pixel assegnati all’oggetto che in realtà non ne fanno parte. I veri positivi (TP) invece saranno i pixel dell’oggetto che effettivamente vengono riconosciuti. Infine, ci sono anche i falsi negativi (FN) ovvero i pixel dell’oggetto che vengono persi. Sulla base di queste tre quantità si possono definire due metriche principali.

C M=\frac{T P}{T P+F N}=\frac{T P}{\text { Total area in GT }}\\
\:\\
C R=\frac{T P}{T P+F P}=\frac{T P}{\text { Total area in BW }}

La completezza (CM) indica quanto la segmentazione ha effettivamente preso dell’oggetto di partenza mentre la correttezza (CR) tiene conto della sovrasegmentazione. Sulla base di questi due si introduce il parametro F-measure (FM) che tiene conto sia della sotto segmentazione che della sovra segmentazione. Tale parametro vive nel range [0; 1] ed è tanto migliore quanto più vicino a 1.

F M=\frac{2 \cdot C M \cdot C R}{C M+C R}\quad  \in \left[0;1\right]

Deep learning

Oggi giorno le reti neurali e l’apprendimento automatico sono molto sviluppate visto anche il recente sviluppo permesso dall’aumento delle potenze di calcolo. Gli ultimi sviluppi che hanno portato all’introduzione sul mercato consumer di computer molto potenti e fotocamere ad alta risoluzione hanno portato all’utilizzo di reti neurali per elaborare le immagini.

L’apprendimento automatico nasce in contrapposizione alla programmazione tradizionale nell’ottica di mettere la macchina nelle condizioni di estrarre autonomamente alcune features di interesse dai dati. A tale scopo nasce quello che prende il nome di deep learning. Diverse tecniche basate su reti artificiali neurali prevedono l’organizzazione su diversi strati ed ognuno di essi calcola dei valori che a loro volta vengono passati al successivo.

Convolutional networks

Dati visuali e più generici dati bidimensionali vengono spesso elaborati con reti neurali convoluzionali, CNN. Sono reti composte da uno o più strati convoluzionali con un’architettura feed-forward. Ovvero le diverse connessioni tra le unità non formano cicli ma le informazioni si muovono su una sola direzione rispetto i nodi di ingresso.

Quindi una CNN non è altro che un algoritmo che prende in ingresso un’immagine e fornendo una certa importanza ad alcuni suoi aspetti. Analizzando diversi aspetti di diverse immagini è in grado di distinguere l’una dalle altre. Le diverse reti CNN hanno alcuni punti in comune nella loro architettutra.

Partendo dal canale di ingresso questo accetterà immagini di una certa dimensione e con un certo numero di canali. Successivamente l’immagine viene passata ad un primo livello convoluzionale dove avviene un primo filtraggio. Questo porta ad estrarre features di alto livello come ad esempio bordi i bruschi cambi di luminosità.

Successivamente ci sono dei layer di polling dove vengono preservate le features che risultano invarianti e posizione e rotazione andando anche a filtrare parte del rumore. I livelli ReLU (rectified Linear Units) hanno l’obiettivo di annullare ed eliminare i valori non utili dopo che i dati sono passati per un livello convoluzionale. Ci sono poi i layer FC (fully connected) che sono quelli che si occupano della classificazione vera e propria.

Transfer learning

È evidente la complessità delle rete neurali e infatti un addestramento accurato richiede un’elevata complessità computazionale. Esiste però una tecnica per adattare un’intelligenza artificiale ad un task diverso da quello per cui è già stata addestrata. L’idea è quella di prende gli strati inferiori di una rete già addestrata e aggiungere nuovi strati finali adattandoli al problema di interesse. Si utilizzando quindi nuovi set di training e algoritmi di ottimizzazione per adattare la classificazione al nuovo problema.

Addestramento

A seguire un’analisi su un’implementazione base all’interno di Matlab di alcune reti neurali per la segmentazione di cellule. Verranno valutate anche le prestazioni di addestramento e di segmentazione. All’interno di Matlab è disponibile il Deep Learning Toolbox che fornisce una piattaforma per la progettazione e l’implementazione di reti neurali profonde con algoritmi e modelli pre addestrati.

Dataset

Il primo passo per addestrare correttamente una rete neurale è quello di selezionare le immagini, le classi e le maschere di identificazione di ciò che si vuole rappresentare. Ovvero si devono considerare le immagini di riferimento e disporre di una maschera che permette di identificare gli oggetti di interesse al suo interno. In questo caso parliamo di un problema binario, cellula e background, quindi si considerano maschere binarie, 0 sul background e 1 sul background.

A livello computazionale in Matlab si riferisce il tutto ad una variabile di datastore. Ovvero un’unica variabile che al suo interno contiene informazioni sul nome dei file, il loro percorso, ecc senza importarli nel workspace e quindi senza occupare memoria. L’intera collegazione verrà poi processata e caricata solo al momento dell’effettivo addestramento.

Training

Il primo passo per poter utilizzare efficacemente una rete neurale preaddestrata è quello di addestrarla correttamente. Addestrare una rete consiste nel fornire un numero adeguato di coppie contenenti l’immagine di test e il risultato della segmentazione. Questo permette alla rete di modificare diversi parametri e coefficienti così da poter svolgere i compiti richiesti.

In particolare, in questo caso si vuole fare una segmentazione binaria, separando il foreground dal background. Quindi l’obiettivo è quello di identificare la cellula e distinguerla dallo sfondo. Nelle prime analisi viene utilizzato un dataset di riferimento diviso in due grandi insiemi. Una prima parte viene utilizzata per l’addestramento e una seconda parte per una successiva analisi della rete.

ResNet50

La prima rete che verrà testata è ResNet50, una rete convoluzionale con 50 layers di profondità. La rete pre addestrata è in grado di classificare oltre 1000 categorie di oggetti. Presenta un totale di 177 layer. Dal layer di input è possibile anche vedere che la rete richiede immagini a 3 canali, rgb, con dimensione 224 x 224. Questo sarà da tenere presente quando viene preparato il dataset di addestramento/test. Per estrarre ulteriori informazioni è possibile andare a leggere la struttura della rete.

net=resnet50;
analyzeNetwork(net)

In alternativa è possibile aprire il Deep Network Analyzer ed esplorare i differenti layer nella più evidente complessità della rete.

Per addestrare la rete il primo passo è considerare il dataset di riferimento e ciò che si vuole ottenere dalla rete. In questo caso si vogliono identificare due classi (background e foreground). In questo caso la rete viene addestrata con 61 immagini inserite all’interno di un ImageDatastore. Ci sono quindi altrettanti ground truth inseriti all’interno di un altro data store, pxds, quest’ultimo però contiene, oltre alle matrici/immagini di GT, anche le informazioni sulla possibilità di avere due classi etichettate con bianco e nero, ovvero rispettivamente 1 e 0.

filePath = matlab.desktop.editor.getActiveFilename;
pathDivided=strsplit(filePath,'\');
newPath=erase(filePath,pathDivided(end));
dataPath=strcat(newPath,'dataset');
addpath(strcat(newPath,'functions')); %set path for functions

Inoltre, è buona norma pre-processare queste immagini affinchè siano compatibili con i layer di input della rete. Il primo passaggio consiste nel generare un dataset aumentato. Ovvero un nuovo ImageDatastore contenente più immagini generate da trasformazioni rigide, rotazioni e crop dell’insieme di partenza. Per farlo si utilizzano i comandi imageDataAugmenter() e il corrispettivo pixelLabelImageDatastore(). Il primo contiene le informazioni sulle operazioni per aumentare il dataset mentre il secondo contiene le informazioni sui dataset di partenza, operazioni di aumento, dimensioni e canali delle immagini, compatibili con quanto richiesto dalla rete.

imds = imageDatastore(strcat(newPath,'dataset\FRAME_TRAIN')); 
pxds = pixelLabelDatastore(strcat(newPath,'dataset\GT_TRAIN'),["N","B"],[0 1]);

Rete convoluzionale

Il passo successivo è quello implementare la rete. Si procede con la creazione di una rete convoluzionale DeepLab v3 direttamente pronta per la segmentazione. Nella sua creazione è necessario specificare la dimensione delle immagini, il numero di classi e la rete di partenza. A questo punto la rete presenterà ulteriori livelli che verrano addestrati.

Tale rete convoluzionale provvederà a classificare i singoli pixel dell’immagine sulla base delle classi con le quali è stata addestrata. Conviene, lavorando con classi binarie dove c’è una predominanza di pixel neri rappresentanti il foreground, normalizzare i coefficienti del layer di segmentazione rispetto la numerosità di ogni classe. Trascurare questa normalizzazione porterebbe ad uno sbilanciamento delle classi e ad un addestramento a favore delle classi predominanti.

numClasses=2; % foreground and background
imageSize=net.Layers(1).InputSize; %read size directly from net
augmenter = imageDataAugmenter('RandRotation',[0 360],'RandXReflection',true,'RandXTranslation',[-10 10],'RandYTranslation',[-10 10]);
pximds = pixelLabelImageDatastore(imds,pxds,'DataAugmentation',augmenter,'OutputSize',imageSize,'ColorPreprocessing','gray2rgb');
lgraph = deeplabv3plusLayers(imageSize, numClasses, "resnet50"); 
% queste blocco serve per bilanciare rispetto alle classi background e
% foreground
tbl = countEachLabel(pximds);
totalNumberOfPixels = sum(tbl.PixelCount);
frequency = tbl.PixelCount / totalNumberOfPixels;
classWeights = 1./frequency;
pxLayer = pixelClassificationLayer('Name','labels','Classes',tbl.Name,'ClassWeights',classWeights);
lgraph = replaceLayer(lgraph,"classification",pxLayer);

Training options

Allora è possibile addestrare la rete. Si utilizza il comando trainNetwork() al quale vengono passati:

• Data store con i dati per l’addestramento
• Layers della rete da addestrare
• Opzioni di addestramento

Quindi è possibile personalizzare alcune opzioni da addestramento. I parametri sicuramente più importanti sono il solutore, il numero di epoche, la dimensione del mini batch, il numero di iterazioni e alcune impostazioni di output.

Il solutore è colui che si occupa di ottimizzare la rete. In Matlab sono implementati diverse varianti di discesa stocastica del gradiente, una tecnica iterativa che si occupa di addestrare la rete. Sono disponibili tre varianti: adam, rmsprop e sdgm, che si comportano in modo diverso e portano a risultati diversi a seconda del problema considerato.

Un’altra variabile di interesse è dal dimensione del batch-size ovvero il sotto insieme sul quale la rete aggiorna i parametri ad ogni iterazione. Ogni ciclo completo invece prende il nome di epoca e a sua volta ne può essere definito il limite superiore, variabile a seconda della dimensione della rete e di quanto è necessario raffinare il transfer learning.

Valutazioni prestazioni

È possibile anche visualizzare l’andamento e il progresso dell’addestramento con l’opzione "Plots","training-progress". Quindi è possibile valutarne l’accuratezza e la perdita dell’addestramento.

options = trainingOptions('sgdm', ...
    'MaxEpochs',30, ...  
    'MiniBatchSize',8, ...
    'Plots','training-progress');

Il comando trainNetwork() lavora principalmente con la GPU e questo permette di velocizzare l’addestramento e scalarne le potenzialità su sistemi multi-GPU. A conferma sono state valutate le prestazioni su due sistemi molto differenti, sia per architettura che per periferiche.

• Mac. Software: MacOS Monterey. Hardware: 3,1 GHz i5, 8 GB DDR3, SSD NVME, Iris + Graphics 650
• Win1. Software: Windows11. Hardware: 2,7 GHz i7-10850H , 32 GB DDR4, SSD NVME, Nvidia Quadro P620
• Win2. Software: Windows11. Hardware: 2,3 GHz i7-11800H, 32 GB DDR4, SSD NVME, Nvidia RTX 3060 Laptop

Vista la grande disparità di velocità e le potenzialità nello sfruttare i driver CUDA su schede Nvidia i risultati che seguono fanno riferimento alla workstation windows Win1.

Test

Una prima analisi può essere fatta, a partire dalle metriche di segmentazione su un nuovo dataset di immagini. Caricandolo, passandolo alla rete per segmentare e valutando il parametro FM è possibile avere una stima pratica di quanto funziona bene la rete.

[net, info]= trainNetwork(pximds,lgraph,options);
accuracy=info.TrainingAccuracy(end);
loss=info.TrainingLoss(end);

for l = 1:length(f_test)
    l
testImage=imread([strcat(dataPath,'/FRAME_TEST_SEG/'),f_test(l).name]);
C_test = semanticseg(testImage,net);
D=C_test=='B';
GTImage=imread([strcat(dataPath,'/GT_TEST/'),gt_train(l).name]);
[TP,FP,FN,CR,CM,FM_test(l)]=evaluation_segmentation(bwareafilt(D,1),GTImage);
imshowpair(testImage,bwareafilt(D,1),'montage');
pause(0.5); drawnow;
clear C_test D testImage;
end

Quindi è possibile andare a vedere come alcuni parametri di training influenzano il risultato finale e le prestazioni della rete.

Variazioni solutore

Un primo fattore di interesse è sicuramente il solutore. Andando a testare, a parità di condizione, i differenti solutori è evidente come si ottiene un’accuratezza maggior con RMSProp mentre una perdita minore con SGDM.

Nonostante i parametri di addestramento siano a favore di RMSProp le metriche di segmentazione, ottenuto da un’effettiva azione della rete, si mostrano a favore di SGDM e RMSProp con valori medi maggiori. Inoltre, SGDM presenta, non solo un valore medio maggiore, anche un valore di picco più elevato ed una deviazione standard più contenuta.

I seguenti test sono effettuati a partire da queste condizioni.

Epoche

Testando invece il numero di epoche è evidente come un valore troppo basso rende la rete poco precisa. Nonostante con un massimo di 60 epoche si ottiene un’accuratezza maggiore, si ottiene un FM ugualmente alto con con solo 20 epoche. Questo a fronte di un numero di iterazioni di 140 contro 420, ovvero una maggior velocità di addestramento.

Dimensione del batch

C’è anche la possibilità di analizzare la dimensione minima del batch che mostra come la scelta iniziale sia la scelta ottima.

Learning rate

Analizzando la velocità iniziale di apprendimento risultano evidenti i valori ottimi. Conviene scegliere quindi il valore 0.01 che porta ad valori di FM più elevati e meno distribuiti.

Augmenter

L’addestramento si basa prevalentemente sulla bontà e grandezza dell’insieme di training. Esistono pertanto alcuni escamotage per ottenere un maggior numero immagini partendo da quelle di partenza e applicando rotazioni, traslazioni e riflessioni. I risultati migliori si hanno grazie alla riflessione, infatti i risultati con solo la riflessione presentano medie di FM molto più alte.

Inoltre, indagando sulla traslazione è possibile osservare come le prestazioni migliori si hanno con traslazioni di +/- 50 pixel.

Considerando tutte le opzioni disponibili, quindi la simultanea rotazione, riflessione e traslazione, si ottengono dei risultati più stabili e comunque elevati. La finestra di traslazione ottima però si sposta a valori più alti per +/- 20 pixel. Tuttavia, dai grafici è evidente come le variazioni siano piccole e la F-measure rimane comunque molto elevato.

Librerie e reti pre addestrate

Si può indagare anche l’evoluzione e la variazione delle metriche di segmentazione andando a sfruttare differenti reti. In particolare l’algortimo Deep Net V3+ permette di utilizzare cinque differenti reti:

• ResNet50
• ResNet18
• InceptionResNetV2
• Xception
• MobileNetV2

Dai risultati è evidente come la scelta migliore sia ResNet50. Le altre reti hanno risultati meno accurati e meno di InceptionResNetV2 che tuttavia presenta un tempo di calcolo 5 volte superiore.

Segmentazione del bordo

Un altro topic di interesse è la ricerca di una regione di bordo di una particolare cellula. Tale ricerca è facilmente estensibile ad un qualsiasi oggetto presente nell’immagine. A tal proposito, e per mostrare la scalabilità di tali algoritmi, viene considerato un dataset di nei per la dermatoscopia. Sebbene possano sembrare immagini molto diverse il contenuto informativo dell’imagine è il medesimo: un oggetto e il suo bordo immersi in un background privo di interesse.

pathImage = strcat(newPath,'dataset\mole\Immagini'); % link for training frame
pathGT=strcat(newPath,'dataset\mole\Segmentazioni');
pathExport=strcat(newPath,'dataset\mole\PROCESSED\');
[imds,pxds]=optimizeDataset(pathImage,pathGT,pathExport,8,imageSize);

Partendo dal dataset vengono preparate le immagini ridimensionate per renderle compatibili con la rete. Viene inoltre preparato il dataset di training facendo in modo di aver 0 sul background, 1 sul bordo e 2 all’interno dell’oggetto.

pxds = pixelLabelDatastore(strcat(pathExport,'GT'),["BACK","BORDER","INSIDE"],[0 1 2]); % link for GT images
imds = imageDatastore(strcat(pathExport,'IMG'));

Addestrando e testando la rete è possibile verificarne il corretto funzionamento su un differente insieme di test.

[area]=testBorderSegmentation(imdsTEST,1,net,'NO');

Trasformazioni affini

In questo caso, non avendo più un problema binario e avendo 3 canali (rgb) cambia leggermente la creazione dell’AugmentedDataStore. Il principio rimane lo stesso del caso precedente ma utilizza funzioni differenti.

dsTrain = combine(imds, pxds);
xTrans = [-10 10];
yTrans = [-10 10];
dsTrain = transform(dsTrain, @(data)augmentImageAndLabel(data,xTrans,yTrans));
lgraph = deeplabv3plusLayers(imageSize, numClasses, "resnet50"); 

Dopo aver combinato il datastore con le immagini e quello con le labels dei singoli pixels va ad utilizzare una funzione di trasformazione che permette di traslare, ridimensionare, ruotare, ecc. le singole immagini.

function data = augmentImageAndLabel(data, xTrans, yTrans)
for i = 1:size(data,1)
    tform = randomAffine2d(...
        'Rotation',[0 360],...
        'XReflection',true,...
        'XTranslation', xTrans, ...
        'YTranslation', yTrans);
    rout = affineOutputView(size(data{i,1}), tform, 'BoundsStyle', 'centerOutput');
    data{i,1} = imwarp(data{i,1}, tform, 'OutputView', rout);
    data{i,2} = imwarp(data{i,2}, tform, 'OutputView', rout);
end
end

Operatori morfologici

Vengono utilizzati anche alcuni operatori morfologici per generare il bordo a partire dal GT che identifica l’intero oggetto. Viene effettuata un’erosione sul GT e poi viene sommato al precedente bordo. Questo porta, passando da operatori logici a numeri interi, ad avere:

• 0 sul background con label 'BACK'
• 1 sul bordo dell’oggetto con label 'BORDER'
• 2 all’interno dell’oggetto con label 'INSIDE'

Per la presentazione finale viene inoltre utilizzato un’operatore di chiusura su ciò che viene identificato come interno dell’oggetto. Quindi l’interno dell’oggetto viene sottratto al bordo così da ripulire ulteriormente l’immagine. Ottimizzando le dimensioni degli elementi con cui agiscono gli operatori morfologici è possibile aumentare la precisione dell’algoritmo.

Conclusioni

La segmentazioni con reti convoluzionali presenta grandi potenziali nel campo medicale. Risulta utile in biologia, una volta identificata la cellula è possibile utilizzare le informazioni in essa contenute per qualsivoglia analisi. Allo stesso tempo è facilmente scalabile all’identificazione di particolari oggetti quali tumori, macchine, fratture, alterazioni, ecc che possono essere presenti all’interno delle diverse immagini medicali.

Codice

Il codice è presente alla repository Git Hub: https://github.com/mastroalex/neural-network

Fonti

• Image Processing Using Deep Learning – Mathworks
• Pretrained Convolutional Neural Networks – Mathworks
• Encoder-Decoder with Atrous Separable Convolution for Semantic Image Segmentation – Liang-Chieh Chen
• Artificial Convolutional Neural Network in Object Detection and Semantic Segmentation for Medical Imaging Analysis – Ruixin Yang
• Deep Semantic Segmentation of Natural and Medical Images: A Review – Saeid Asgari Taghanaki
• Deep Neural Architectures for Medical Image Semantic Segmentation: Review – Muhammad Zubair Khan