Vai al contenuto principale

Tutorial: Super-risoluzione di immagini con framework Multi-Decoder

In questo tutorial, implementerai un paper di imaging medico con deep learning in Python usando Keras.
Aggiornato 29 set 2026  · 15 min leggi

Scopri con l'IA

ChatGPTClaudePerplexity

Nella parte dedicata allo script di training, lavorerai sul problema della super-risoluzione di immagini usando una nuova architettura di deep learning. Il compito consiste nel mappare in modo non lineare un’immagine RM cerebrale a basso campo 3 Tesla in un’immagine RM cerebrale ad alto campo 7 Tesla. In questa parte dedicata allo script di testing, userai i pesi addestrati nella Parte 1 per fare previsioni su dati mai visti in addestramento. Imparerai anche come salvare infine le immagini 2D come un volume combinato.

Il tutorial è diviso in due parti: la prima ti guiderà attraverso il processo di training e la seconda coprirà il processo di testing.

Nota: Probabilmente ti interesserà leggere il paper: puoi trovare l’articolo qui.

Training Script

In sintesi, nella Parte 1 del tutorial affronterai i seguenti temi:

  • Inizierai importando i moduli necessari per addestrare il tuo modello di deep learning,
  • Poi una breve introduzione al dataset MRI 3T e 7T,
  • Quindi definirai gli initializer e caricherai i dataset 3T e 7T; durante il caricamento ridimensionerai anche le immagini al volo,
  • Successivamente preprocesserai i dati caricati: convertirai le liste di train e test in matrici numpy, convertirai il tipo delle matrici in float32, riscalerai le matrici con una strategia max-min, rimodellerai gli array e infine dividerai i dati in 80% training e 20% validation,
  • Poi creerai l’architettura 1-Encoder-3-Decoder: con connessioni di merge e multi-decoder,
  • Successivamente definirai la loss function, tre modelli diversi e li compilerai,
  • Infine, è il momento di addestrare il tuo modello con merge e multi-decoder, testarlo sui dati di validazione e calcolare i risultati quantitativi.

Dipendenze dei moduli Python

Prima di seguire questo tutorial, assicurati di avere esattamente le stesse versioni dei moduli indicate di seguito:

  • Keras==2.0.4
  • tensorflow==1.8.0
  • scipy==0.19.0
  • numpy==1.14.5
  • Pillow==4.1.1
  • nibabel==2.1.0
  • scikit_learn==0.18.1

Nota: Prima di iniziare, considera che il modello verrà addestrato su un sistema con GPU Nvidia 1080 Ti, processore Xeon e5 GeForce e 32 GB di RAM. Se usi Jupyter Notebook, dovrai aggiungere tre righe di codice in cui specifichi l’ordine dei device CUDA e i dispositivi CUDA visibili usando il modulo os.

Nel codice seguente imposti variabili d’ambiente nel notebook con os.environ. È consigliabile farlo prima di inizializzare Keras per limitare TensorFlow (backend di Keras) a usare la prima GPU. Se la macchina su cui addestri ha una GPU su 0, assicurati di usare 0 invece di 1. Puoi verificarlo eseguendo un semplice comando nel terminale: ad esempio, nvidia-smi

import os
os.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"
os.environ["CUDA_VISIBLE_DEVICES"]="0" #model will be trained on GPU

Importing the modules

Per prima cosa, importa tutti i moduli necessari come tensorflow, numpy e soprattutto keras e le funzioni o layer richiesti come Input, Conv2D, MaxPooling2D ecc., poiché userai questi framework per addestrare il modello!
Per leggere immagini in formato NIfTI, devi anche importare un modulo chiamato nibabel.

import os
import numpy as np
import math
import tensorflow as tf
import nibabel as nib
import numpy as np
from keras.layers import Input,Dense,merge,Reshape,Conv2D,MaxPooling2D,UpSampling2D
from keras.layers.normalization import BatchNormalization
from keras.models import Model,Sequential
from keras.callbacks import ModelCheckpoint
from keras.optimizers import RMSprop
from keras import backend as K
import scipy.misc
from sklearn.utils import shuffle
from sklearn.cross_validation import train_test_split
import matplotlib.pyplot as plt
from keras.models import model_from_json

Using TensorFlow backend.
/usr/local/lib/python3.5/dist-packages/sklearn/cross_validation.py:41: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also, note that the interface of the new CV iterators is different from that of this module. This module will be removed in 0.20.
  "This module will be removed in 0.20.", DeprecationWarning)

Understanding the Brain MRI 3T and 7T Dataset

Il dataset brain MRI 3T e 7T è composto da volumi 3D, ciascun volume ha in totale 207 slice/immagini RM cerebrali acquisite a diverse sezioni del cervello. Ogni slice ha dimensione 173 × 173. Le immagini sono in scala di grigi a canale singolo. Ci sono in totale 39 soggetti, ciascuno contenente la scansione RM di un paziente. Il formato delle immagini non è jpeg, png, ecc., ma il formato NIfTI. Nella sezione successiva vedrai come leggere le immagini in formato NIfTI.

Il dataset consiste in immagini RM di tipo T1; le sequenze T1 sono tradizionalmente considerate adatte per la valutazione delle strutture anatomiche. Il dataset su cui lavorerai oggi consiste in RM cerebrali 3T e 7T.

Il dataset è pubblico ed è disponibile per il download a questa fonte.

28 soggetti sono usati per il training mentre i restanti 11 soggetti verranno usati per il testing.

Defining the Initializers

Definiamo innanzitutto le dimensioni dei dati. Ridimensionerai le immagini da 173×173 a 176×176 nella parte di lettura dei dati. Qui definirai anche la directory dei dati, la batch size per l’addestramento, il numero di canali, un layer Input(), le matrici di train e test come liste e infine, per il riscalamento, caricherai il file di testo che contiene i valori minimo e massimo del dataset MRI.

Nota: Per il riscalamento, puoi anche calcolare massimo e minimo del dataset, dato che potresti non avere il file maxANDmin.txt.

x,y = 173,173
full_z = 207
resizeTo=176
batch_size = 32
inChannel = outChannel = 1
input_shape=(x,y,inChannel)
input_img = Input(shape = (resizeTo, resizeTo, inChannel))
inp = "ground3T"
out = "ground7T"     
train_matrix = []
test_matrix = []
min_max = np.loadtxt('maxANDmin.txt')

Loading the Data

Ora carichi i dati RM usando la libreria nibabel e ridimensioni le immagini da 173 × 173 a 176 × 176, aggiungendo zeri in padding nelle dimensioni x e y.

Nota che quando carichi un volume in formato NIfTI, Nibabel non carica subito l’array dell’immagine: attende finché non richiedi i dati dell’array. Il modo standard è chiamare il metodo get_data().

Poiché ti servono slice 2D invece del 3D, userai le liste train e test inizializzate prima; ogni volta che leggi un volume, iteri su tutte le 207 slice del volume 3D e accodi ciascuna slice una per una in una lista.

folder = os.listdir(inp)

Carichiamo prima i dati 3T.

for f in folder:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    a = nib.load(inp + f)
    a = a.get_data()
    temp[3:,:,3:] = a
    a = temp
    for j in range(full_z):
        train_matrix.append(a[:,j,:])

Poi carica i dati 7T, e usa la stessa variabile folder anche per i 7T poiché il numero di volumi RM 3T e 7T è uguale.

for f in folder:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    b = nib.load(out + f)
    b = b.get_data()
    temp[3:,:,3:] = b
    b = temp
    for j in range(full_z):
        test_matrix.append(b[:,j,:])

Data Preprocessing

Dato che le matrici di train e test sono liste, userai numpy per convertirle in array numpy.

Inoltre convertirai il type dell’array numpy in float32 e riscalerai sia l’input sia il ground truth.

train_matrix = np.asarray(train_matrix)
train_matrix = train_matrix.astype('float32')
m = min_max[0]
mi = min_max[1]
train_matrix = (train_matrix - mi) / (m - mi)

test_matrix = np.asarray(test_matrix)
test_matrix = test_matrix.astype('float32')
test_matrix = (test_matrix - mi) / (m - mi)

Stampiamo velocemente la shape delle matrici train_matrix (3T) e test_matrix (7T). Dovrebbero avere in totale 28 × 207 = 5796 immagini, ciascuna con dimensione 176 × 176.

train_matrix.shape

(5796, 176, 176)

test_matrix.shape

(5796, 176, 176)

Successivamente creerai due nuove variabili augmented_images (3T/input) e Haugmented_images (7T/ground truth) con la stessa shape di train e test. Saranno matrici 4D, in cui la prima dimensione è il numero totale di immagini, la seconda e la terza sono le dimensioni di ogni immagine e l’ultima è il numero di canali, che in questo caso è uno.

augmented_images=np.zeros(shape=[(train_matrix.shape[0]),(train_matrix.shape[1]),(train_matrix.shape[2]),(1)])
Haugmented_images=np.zeros(shape=[(train_matrix.shape[0]),(train_matrix.shape[1]),(train_matrix.shape[2]),(1)])

Poi itera su tutte le immagini una per una; ogni volta rimodella le matrici di train e test a 176 × 176 e le accoda rispettivamente in augmented_images (3T/input) e Haugmented_images (7T/ground truth).

for i in range(train_matrix.shape[0]):
    augmented_images[i,:,:,0] = train_matrix[i,:,:].reshape(resizeTo,resizeTo)
    Haugmented_images[i,:,:,0] = test_matrix[i,:,:].reshape(resizeTo,resizeTo)

Dopo tutto ciò, è importante partizionare i dati. Per favorire la capacità di generalizzazione del modello, dividi i dati in due parti: training e validation. Addestrerai il modello sull’80% dei dati e lo validerai sul restante 20%.

Questo aiuta anche a ridurre le probabilità di overfitting, poiché validerai il modello su dati non visti in fase di training.

Puoi usare il modulo train_test_split di scikit-learn che hai importato all’inizio per dividere correttamente i dati:

data,Label = shuffle(augmented_images,Haugmented_images, random_state=2)
X_train, X_test, y_train, y_test = train_test_split(data, Label, test_size=0.2, random_state=2)

X_test = np.array(X_test)
y_test = np.array(y_test)
X_test = X_test.astype('float32')
y_test = y_test.astype('float32')

X_train = np.array(X_train)
y_train = np.array(y_train)
X_train = X_train.astype('float32')
y_train = y_train.astype('float32')

Il modello: 1-Encoder-3-Decoders!

model
Figura: Backpropagation dell’autoencoder selettivo
Immagine tratta da questo paper.

Merge Connections

Ora definisci l’architettura proposta con blocchi di layer di filtri consecutivi seguiti da un layer di max pooling nella sezione di encoder, come mostrato nella cella successiva. Per ricostruire la dimensione originale dell’immagine in uscita, in ciascun blocco dei decoder è introdotto un layer di upsampling. Durante l’upsampling, possono emergere artefatti dovuti ai dettagli mancanti nell’input downsampled del decoder. Pertanto, concateni l’input del decoder con la sua versione upscalata dall’encoder, in modo da fornire la natura dei dettagli upscalati per una migliore ricostruzione durante l’upsampling nel decoder. L’aggiunta delle merge connections produce infatti un miglioramento significativo del PSNR (dell’ordine di 5 dB). Questa impostazione dell’architettura è ispirata a questo paper.

Multi Decoders

L’approccio proposto impiega un singolo encoder e un’architettura con più decoder con input a canale singolo. In ciascun blocco dell’encoder e in tutti e tre i decoder si usano tre layer convoluzionali, seguiti da un layer di batch normalization per mantenere la stabilità numerica.

Il primo blocco convoluzionale dell’encoder ha 32 filtri, e il numero di filtri raddoppia dopo ogni blocco convoluzionale. In tutti i decoder, il primo blocco convoluzionale ha 256 filtri, e il numero di filtri si dimezza dopo ogni blocco. Userai un filtro 3×3 in tutti i blocchi convoluzionali.

La Rectified Linear Unit (ReLU) sarà usata come funzione di attivazione in tutti i layer tranne l’ultimo. Poiché i dati sono normalizzati tra 0 e 1, nell’ultimo layer si usa una Sigmoid.

È noto che i dettagli locali dell’immagine a varie scale svolgono un ruolo significativo nella ricostruzione. L’architettura proposta considera immagini a scale diverse usando layer gerarchici per il downsampling (max pooling) e l’upsampling (entrambi con fattore 2) rispettivamente in encoder e decoder. La rappresentazione codificata ottenuta dopo tre operazioni di downsampling porta i dati da un input ad alta dimensione a una rappresentazione in spazio latente.

Encoder

def encoder(input_img):
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(input_img)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(pool3)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv4)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='sigmoid', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    return conv5,conv4,conv3,conv2,conv1

Decoder 1

def decoder_1(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded_1 = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded_1

Decoder 2

def decoder_2(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded_2 = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded_2

Decoder 3

def decoder_3(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded_3 = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded_3

Nelle 4 celle sopra hai definito quattro funzioni: una per l’encoder e tre per i decoder. Trattandosi di funzioni, potresti definire una funzione decoder() e chiamarla tre volte; ma per migliore comprensione, le definiremo tre volte e per evitare problemi dovuti alla casualità dei tre decoder.

Loss Function

Ora userai un mean square error in cui escluderai i valori (pixel) sia del ground truth y_t sia di y_p uguali a zero.

def root_mean_sq_GxGy(y_t, y_p):
    a1=1
    where = tf.not_equal(y_t, 0)
    a_t=tf.boolean_mask(y_t,where,name='boolean_mask')
    a_p=tf.boolean_mask(y_p,where,name='boolean_mask')
    return a1*(K.sqrt(K.mean((K.square(a_t-a_p)))))

Model Definition and Compilation

Per prima cosa, chiami la funzione encoder passandole l’input. Poiché usi le merge connection nell’architettura, la funzione encoder restituisce l’output di cinque layer di convolution, che poi mergi separatamente con l’output di tutti e tre i decoder.

conv5,conv4,conv3,conv2,conv1 = encoder(input_img)

autoencoder_1 = Model(input_img, decoder_1(conv5,conv4,conv3,conv2,conv1))
autoencoder_1.compile(loss=root_mean_sq_GxGy, optimizer = RMSprop())

autoencoder_2 = Model(input_img, decoder_2(conv5,conv4,conv3,conv2,conv1))
autoencoder_2.compile(loss=root_mean_sq_GxGy, optimizer = RMSprop())

autoencoder_3 = Model(input_img, decoder_3(conv5,conv4,conv3,conv2,conv1))
autoencoder_3.compile(loss=root_mean_sq_GxGy, optimizer = RMSprop())

/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:2: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.

/usr/local/lib/python3.5/dist-packages/keras/legacy/layers.py:460: UserWarning: The `Merge` layer is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.
  name=name)
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:10: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.
  # Remove the CWD from sys.path while we load stuff.
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:18: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:26: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.


WARNING:tensorflow:From /usr/local/lib/python3.5/dist-packages/keras/backend/tensorflow_backend.py:1257: calling reduce_mean (from tensorflow.python.ops.math_ops) with keep_dims is deprecated and will be removed in a future version.
Instructions for updating:
keep_dims is deprecated, use keepdims instead

Train the Model

Salverai i pesi solo quando il peak signal-to-noise ratio sui dati di validazione migliora. A tal fine definirai una lista psnr_gray_channel in cui accodare un valore predefinito pari a 1. Inserisci 1 perché devi inizializzare la lista con un numero qualsiasi, e il PSNR sarà molto più alto anche nelle fasi iniziali di training: fungerà solo da numero fittizio.

Inizializzerai il learning rate a 1e-3 e userai una strategia di decay in cui lo riduci del 10% rispetto al valore corrente ogni 20 epoche.

psnr_gray_channel = []
psnr_gray_channel.append(1)
learning_rate = 0.001
j=0

Nota: Il prossimo blocco di codice andrebbe eseguito in un’unica cella, ma per una migliore comprensione dell’intero processo di training verrà suddiviso in più celle.

Il modello è addestrato per 500 epoche. Il learning rate iniziale è 1e-3 come definito in precedenza. Salverai anche i valori di PSNR e MSE dopo ogni epoca tra 7T e il predetto. Userai la funzione K.set_value per modificare il learning rate di tutti e tre i modelli del 10%, ma ogni 20 epoche.

for jj in range(500):
    myfile_valid_psnr_7T = open('../1_encoder_3_decoders_complete_slices_single_channel/validation_psnr7T_1encoder_3decoders.txt', 'a')
    myfile_valid_mse_7T = open('../1_encoder_3_decoders_complete_slices_single_channel/validation_mse7T_1encoder_3decoders.txt', 'a')

    K.set_value(autoencoder_1.optimizer.lr, learning_rate)
    K.set_value(autoencoder_2.optimizer.lr, learning_rate)
    K.set_value(autoencoder_3.optimizer.lr, learning_rate)

Successivamente rimescoli le immagini 3T di input e le 7T di ground truth per evitare l’overfitting, dato che non rimescolare costringerebbe il modello a vedere i campioni nello stesso ordine a ogni epoca. Poi calcoli il numero di batch che verranno creati in base alla batch_size definita prima. Infine inizi a iterare su num_batches.

train_X,train_Y = shuffle(X_train,y_train)
print ("Epoch is: %d\n" % j)
print ("Number of batches: %d\n" % int(len(train_X)/batch_size))
num_batches = int(len(train_X)/batch_size)
for batch in range(num_batches):

Oltre a memorizzare i valori PSNR, salverai anche le loss di tutti e tre gli autoencoder e dei tre decoder rispettivamente.

myfile_ae1_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/ae1_train_loss_1encoder_3decoders.txt', 'a')
myfile_ae2_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/ae2_train_loss_1encoder_3decoders.txt', 'a')
myfile_ae3_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/ae3_train_loss_1encoder_3decoders.txt', 'a')
myfile_dec1_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/dec1_train_loss_1encoder_3decoders.txt', 'a')
myfile_dec2_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/dec2_train_loss_1encoder_3decoders.txt', 'a')
myfile_dec3_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/dec3_train_loss_1encoder_3decoders.txt', 'a')

Dato che in ogni batch vuoi che il modello veda i successivi 32 campioni (batch_size), la prossima cella si occupa di questo.

batch_train_X = train_X[batch*batch_size:min((batch+1)*batch_size,len(train_X)),:]
batch_train_Y = train_Y[batch*batch_size:min((batch+1)*batch_size,len(train_Y)),:]

Per far funzionare la nostra strategia del Minimum autoencoder, dopo ogni epoca testi sui dati di training usando la funzione Keras test_on_batch, che ti restituirà tre loss rispettivamente, e infine le stampi.

loss_1 = autoencoder_1.test_on_batch(batch_train_X,batch_train_Y)
loss_2 = autoencoder_2.test_on_batch(batch_train_X,batch_train_Y)
loss_3 = autoencoder_3.test_on_batch(batch_train_X,batch_train_Y)
print ('epoch_num: %d batch_num: %d Test_loss_1: %f\n' % (j,batch,loss_1))
print ('epoch_num: %d batch_num: %d Test_loss_2: %f\n' % (j,batch,loss_2))
print ('epoch_num: %d batch_num: %d Test_loss_3: %f\n' % (j,batch,loss_3))

Nel tuo network ci possono essere sei condizioni possibili:
loss_1: loss dell’Autoencoder 1
loss_2: loss dell’Autoencoder 2
loss_3: loss dell’Autoencoder 3

    • loss_1 può essere maggiore di loss_2 e loss_3. Se è vero, allora addestri usando solo l’Autoencoder 1. Imposti la parte di Encoder di Autoencoder 2 e Autoencoder 3 a False e addestri solo i loro decoder, quindi scrivi tutte le loss nei file di testo.
    • loss_2 può essere maggiore di loss_1 e loss_3. Se è vero, allora addestri usando solo l’Autoencoder 2. Imposti la parte di Encoder di Autoencoder 1 e Autoencoder 3 a False e addestri solo i loro decoder, quindi scrivi tutte le loss nei file.
    • loss_3 può essere maggiore di loss_1 e loss_2. Se è vero, allora addestri usando solo l’Autoencoder 3. Imposti la parte di Encoder di Autoencoder 1 e Autoencoder 2 a False e addestri solo i loro decoder, quindi scrivi tutte le loss nei file.
    • loss_1 può essere uguale a loss_2. Se è vero, allora addestri usando uno tra Autoencoder 1 o 2. Imposti la parte di Encoder di Autoencoder 3 a False insieme a quella dell’autoencoder non scelto tra 1 e 2, e addestri solo i loro decoder, quindi scrivi tutte le loss nei file.
    • loss_2 può essere uguale a loss_3. Se è vero, allora addestri usando uno tra Autoencoder 2 o 3. Imposti la parte di Encoder di Autoencoder 1 a False insieme a quella dell’autoencoder non scelto tra 2 e 3, e addestri solo i loro decoder, quindi scrivi tutte le loss nei file.
    • loss_3 può essere uguale a loss_1. Se è vero, allora addestri usando uno tra Autoencoder 3 o 1. Imposti la parte di Encoder di Autoencoder 2 a False insieme a quella dell’autoencoder non scelto tra 3 e 1, e addestri solo i loro decoder, quindi scrivi tutte le loss nei file.
model
Figura: Architettura del modello
Immagine tratta da questo paper.
if loss_1 < loss_2 and loss_1 < loss_3:
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae1_loss.write("%f \n" % (train_1))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_1: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_1: %f\n' % (j,batch,train_1))
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec2_loss.write("%f \n" % (train_2))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
    print ('epoch_num: %d batch_num: %d Decoder_loss_3: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_loss_3: %f\n' % (j,batch,train_3))
elif loss_2 < loss_1 and loss_2 < loss_3:
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae2_loss.write("%f \n" % (train_2))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_2: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_2: %f\n' % (j,batch,train_2))
    for layer in autoencoder_1.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    #autoencoder_1.summary()
    #autoencoder_3.summary()
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec1_loss.write("%f \n" % (train_1))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3: %f\n' % (j,batch,train_3))
elif loss_3 < loss_1 and loss_3 < loss_2:
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_3: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_3: %f\n' % (j,batch,train_3))
    for layer in autoencoder_1.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    #autoencoder_1.summary()
    #autoencoder_2.summary()
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec1_loss.write("%f \n" %(train_1))
    myfile_dec2_loss.write("%f \n" % (train_2))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
elif loss_1 == loss_2:
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae1_loss.write("%f \n" % (train_1))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec2_loss.write("%f \n" % (train_2))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
elif loss_2 == loss_3:
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae2_loss.write("%f \n" % (train_2))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    for layer in autoencoder_1.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False

    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec1_loss.write("%f \n" % (train_1))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
elif loss_3 == loss_1:
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae1_loss.write("%f \n" % (train_1))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec2_loss.write("%f \n" % (train_2))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))


    myfile_ae1_loss.close()
    myfile_ae2_loss.close()
    myfile_ae3_loss.close()
    myfile_dec1_loss.close()
    myfile_dec2_loss.close()
    myfile_dec3_loss.close()

Questo è un passaggio essenziale: avendo impostato alcuni layer a False nelle condizioni di cui sopra, devi reimpostarli a True, così che tutti i layer degli autoencoder possano essere usati dalla funzione test_on_batch e non rimangano False per tutto il training.

for layer in autoencoder_1.layers[:34]:
            layer.trainable = True
for layer in autoencoder_2.layers[:34]:
            layer.trainable = True
for layer in autoencoder_3.layers[:34]:
            layer.trainable = True

Salverai i pesi in tre modi diversi; sotto ne vedi due. Uno: salvi i pesi ogni 100 epoche. Due: salvi i pesi dopo ogni epoca, sovrascrivendoli.

if jj % 100 ==0:
            autoencoder_1.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE1_" + str(jj)+".h5")
            autoencoder_2.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE2_" + str(jj)+".h5")
            autoencoder_3.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE3_" + str(jj)+".h5")


autoencoder_1.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE1.h5")
autoencoder_2.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE2.h5")
autoencoder_3.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE3.h5")

Testing on Validation Data

Per prima cosa rimescoli i dati di validazione e poi usi le stesse 6 possibili condizioni descritte sopra. Qualunque condizione sia vera, testi il modello usando quell’autoencoder specifico.

Quindi calcoli due metriche, MSE e PSNR, tra decoded_imgs (predetto) e ground truth. Infine le salvi nei file di testo.

X_test,y_test = shuffle(X_test,y_test)
if loss_1 < loss_2 and loss_1 < loss_3:
    decoded_imgs = autoencoder_1.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)


    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))

    #print (check)
elif loss_2 < loss_1 and loss_2 < loss_3:
    decoded_imgs = autoencoder_2.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))

    #print (check)
elif loss_3 < loss_2 and loss_3 < loss_1:
    decoded_imgs = autoencoder_3.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))
    #print (check)

elif loss_1 == loss_2:
    decoded_imgs = autoencoder_1.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))


elif loss_2 == loss_3:
    decoded_imgs = autoencoder_2.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))


elif loss_3 == loss_1:
    decoded_imgs = autoencoder_3.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))

Qui salvi i pesi solo quando il PSNR tra il 7T predetto e il ground truth (7T) è massimo rispetto ai precedenti valori di PSNR salvati nella lista psnr_gray_channel.

if max(psnr_gray_channel) < psnr_7T:
            autoencoder_1.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE1_" + str(jj)+".h5")
            autoencoder_2.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE2_" + str(jj)+".h5")
            autoencoder_3.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE3_" + str(jj)+".h5")

    psnr_gray_channel.append(psnr_7T)

Saving Input, Ground Truth and Decoded: Quantitative Results

Definirai una matrice numpy temp di dimensioni 176 × 528, dato che salverai 3 immagini in una riga, ciascuna 176 × 176. Salverai una delle immagini tra tutte quelle di validazione 3T, 7T e predetta, e moltiplicherai la matrice per 255 poiché le immagini erano state scalate tra 0 e 1.

Infine, usando la funzione di scipy, salverai l’immagine.

temp = np.zeros([resizeTo,resizeTo*3])
temp[:resizeTo,:resizeTo] = X_test[0,:,:,0]
temp[:resizeTo,resizeTo:resizeTo*2] = y_test[0,:,:,0]
temp[:resizeTo,2*resizeTo:] = decoded_imgs[0,:,:,0]
temp = temp*255
scipy.misc.imsave('../Results/1_encoder_3_decoders_complete_slices_single_channel/' + str(j) + '.jpg', temp)
j +=1

Chiudiamo i file di PSNR e MSE alla fine.

myfile_valid_psnr_7T.close()
myfile_valid_mse_7T.close()

Infine, alla fine, riduci il learning rate del 10% del valore corrente ogni 20 epoche.

if jj % 20 ==0:
        learning_rate = learning_rate - learning_rate * 0.10

Testing Script

Nella parte 2 di questo tutorial affronterai i seguenti temi:

  • Inizierai importando i moduli necessari per addestrare il tuo modello di deep learning,
  • Poi una breve introduzione al dataset MRI 3T e 7T,
  • Quindi definirai gli initializer e caricherai i dataset di test 3T e 7T; durante il caricamento ridimensionerai anche le immagini al volo,
  • Successivamente preprocesserai i dati caricati: convertirai le liste di train e test in matrici numpy, convertirai il tipo delle matrici in float32, riscalerai le matrici con strategia max-min, rimodellerai gli array e infine dividerai i dati in 80% training e 20% validation,
  • Poi creerai l’architettura 1-Encoder-3-Decoder: con connessioni di merge e multi-decoder,
  • Successivamente definirai: loss function, tre modelli diversi e infine caricherai i pesi addestrati,
  • Infine, è il momento di predire con il tuo modello con merge e multi-decoder sui dati non visti e salvare sia i risultati quantitativi sia quelli qualitativi. Imparerai anche a salvare immagini 2D come un volume combinato usando nibabel.

Dipendenze dei moduli Python

Prima di seguire questo tutorial, assicurati di avere esattamente le stesse versioni dei moduli indicate di seguito:

  • Keras==2.0.4
  • tensorflow==1.8.0
  • scipy==0.19.0
  • numpy==1.14.5
  • Pillow==4.1.1
  • nibabel==2.1.0
  • scikit_learn==0.18.1

Nota: Prima di iniziare, considera che il modello verrà addestrato su un sistema con GPU Nvidia 1080 Ti, processore Xeon e5 GeForce e 32 GB di RAM. Se usi Jupyter Notebook, dovrai aggiungere tre righe di codice in cui specifichi l’ordine dei device CUDA e i dispositivi CUDA visibili usando il modulo os.

Nel codice seguente, imposti sostanzialmente variabili d’ambiente nel notebook con os.environ. È buona pratica farlo prima di inizializzare Keras per limitare TensorFlow (backend di Keras) a usare la prima GPU. Se la macchina su cui addestri ha una GPU su 0, assicurati di usare 0 invece di 1. Puoi verificarlo eseguendo un semplice comando nel terminale: ad esempio, nvidia-smi

import os
os.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"
os.environ["CUDA_VISIBLE_DEVICES"]="0" #model will be trained on GPU

Importing the modules

Per prima cosa, importa tutti i moduli necessari come tensorflow, numpy e soprattutto keras e le funzioni o layer richiesti come Input, Conv2D, MaxPooling2D, ecc., poiché userai questi framework per addestrare il modello!
Per leggere le immagini in formato NIfTI, devi anche importare il modulo nibabel.

import os
from keras.layers import Input,Dense,Flatten,Dropout,merge,Reshape,Conv2D,MaxPooling2D,UpSampling2D
from keras.layers.normalization import BatchNormalization
from keras.models import Model,Sequential
from keras.callbacks import ModelCheckpoint
from keras.optimizers import Adadelta, RMSprop,SGD,Adam
from keras import regularizers
from keras import backend as K
import numpy as np
import scipy.misc
import numpy.random as rng
from sklearn.utils import shuffle
import nibabel as nib
from sklearn.cross_validation import train_test_split
import math

Understanding the Brain MRI 3T and 7T Dataset

Il dataset brain MRI 3T e 7T è composto da volumi 3D, ciascun volume ha in totale 207 slice/immagini RM cerebrali acquisite a diverse sezioni del cervello. Ogni slice ha dimensione 173 × 173. Le immagini sono in scala di grigi a canale singolo. Ci sono in totale 39 soggetti, ciascuno contenente la scansione RM di un paziente. Il formato delle immagini non è jpeg, png, ecc., ma il formato NIfTI. Nella sezione successiva vedrai come leggere le immagini in formato NIfTI.

Il dataset consiste in immagini RM di tipo T1; le sequenze T1 sono tradizionalmente considerate adatte per la valutazione delle strutture anatomiche. Il dataset su cui lavorerai oggi consiste in RM cerebrali 3T e 7T.

Il dataset è pubblico ed è disponibile per il download a questa fonte.

Defining the Initializers

Definiamo innanzitutto le dimensioni dei dati. Ridimensionerai le immagini da 173×173 a 176×176 nella parte di lettura dei dati. Qui definiamo anche la directory dei dati, la batch size usata per l’addestramento, il numero di canali, un layer Input(), le matrici di train e test come liste e infine, per il riscalamento, carichiamo il file di testo che contiene i valori minimo e massimo del dataset MRI.

x,y = 173,173
full_z = 207
resizeTo=176
inChannel = outChannel = 1
input_shape=(x,y,inChannel)
input_img = Input(shape = (resizeTo, resizeTo, inChannel))

train_matrix = []
test_matrix = []
ff = os.listdir("../test_crossval1")
save = "../Result_nii_crossval1"
folder_ground = os.listdir("../test_g_crossval1")
ToPredict_images=[]
predict_matrix=[]

ground_images=[]
ground_matrix=[]
min_max=np.loadtxt('../maxANDmin.txt')

Load the Testing Volumes

Ora carichiamo i dati RM usando la libreria nibabel e ridimensioniamo le immagini da 173 × 173 a 176 × 176 aggiungendo zeri in padding nelle dimensioni x e y.

Nota che quando carichi un volume in formato NIfTI, Nibabel non carica subito l’array dell’immagine: attende finché non richiedi i dati dell’array. Il modo standard è chiamare il metodo get_data().

Poiché vuoi le slice 2D invece del 3D, userai le liste train e test inizializzate prima; ogni volta che leggi un volume, itererai su tutte le 207 slice del volume 3D e accoderai ciascuna slice una per una in una lista.

for f in ff:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    a = nib.load("../test_crossval1" + f)
    affine = a.affine
    a = a.get_data()
    temp[3:,:,3:] = a
    a = temp
    for j in range(full_z):
        predict_matrix.append(a[:,j,:])

for f in ff:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    a = nib.load("../test_g_crossval1" + f)
    affine = a.affine
    a = a.get_data()
    temp[3:,:,3:] = a
    a = temp
    for j in range(full_z):
        ground_matrix.append(a[:,j,:])

Data Preprocessing

Dato che le matrici di test 3T e 7T sono liste, userai numpy per convertirle in array numpy.

Inoltre convertirai il type dell’array numpy in float32 e riscalerai sia l’input sia il ground truth.

ToPredict_images = np.asarray(predict_matrix)
ToPredict_images = ToPredict_images.astype('float32')
mx = min_max[0]
mn = min_max[1]
ToPredict_images[:,:,:,0] = (ToPredict_images[:,:,:,0] - mn ) / (mx - mn)

ground_images = np.asarray(ground_matrix)
ground_images = ground_images.astype('float32')
ground_images[:,:,:,0] = (ground_images[:,:,:,0] - mn ) / (mx - mn)

Poi creerai due nuove variabili ToPredict_images (3T Test/input) e ground_images (7T Test/ground truth) con la stessa shape delle matrici di train e test. Saranno matrici 4D, in cui la prima dimensione è il numero totale di immagini, la seconda e la terza sono le dimensioni di ogni immagine e l’ultima è il numero di canali, che in questo caso è uno.

ToPredict_images=np.zeros(shape=[(ToPredict_images.shape[0]),(ToPredict_images.shape[1]),(ToPredict_images.shape[2]),(1)])
ground_images=np.zeros(shape=[(ground_images.shape[0]),(ground_images.shape[1]),(ground_images.shape[2]),(1)])

Quindi itererai su tutte le immagini una per una; ogni volta rimodellerai le matrici di train e test a 176 × 176 e le accoderai rispettivamente in ToPredict_images (3T Test/input) e ground_images (7T Test/ground truth).

for i in range(ToPredict_images.shape[0]):
    ToPredict_images[i,:,:,0] = ToPredict_images[i,:,:].reshape(resizeTo,resizeTo)

for i in range(ground_images.shape[0]):
    ground_images[i,:,:,0] = ground_images[i,:,:].reshape(resizeTo,resizeTo)

The Model!

model
Figura: Architettura del modello
Immagine tratta da questo paper.

Encoder

def encoder(input_img):
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(input_img)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(pool3)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv4)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='sigmoid', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    return conv5,conv4,conv3,conv2,conv1

Decoder

def decoder(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded

Loss Function

Ora userai un mean square error in cui escluderai i valori (pixel) sia del ground truth y_t sia di y_p uguali a zero.

def root_mean_sq_GxGy(y_t, y_p):
    a1=1
    zero = tf.constant(0, dtype=tf.float32)
    where = tf.not_equal(y_t, zero)
    a_t=tf.boolean_mask(y_t,where,name='boolean_mask')
    a_p=tf.boolean_mask(y_p,where,name='boolean_mask')
    return a1*(K.sqrt(K.mean((K.square(a_t-a_p)))))

Model Definition and Load Weights in all three Autoencoders!

Per prima cosa, richiami la funzione encoder passandole l’input. Poiché usi la merge connection nell’architettura, la funzione encoder restituisce l’output di cinque layer di convolution, che poi mergi separatamente con l’output di tutti e tre i decoder.

conv5,conv4,conv3,conv2,conv1 = encoder(input_img)

Ora creiamo tre modelli diversi e carichiamo in essi i tre pesi addestrati.

autoencoder_1 = Model(input_img, decoder(conv5,conv4,conv3,conv2,conv1))

autoencoder_1.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE1.h5")

autoencoder_2 = Model(input_img, decoder(conv5,conv4,conv3,conv2,conv1))

autoencoder_2.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE2.h5")

autoencoder_3 = Model(input_img, decoder(conv5,conv4,conv3,conv2,conv1))

autoencoder_3.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE3.h5")

Model Prediction on Test Volumes: Quantitative and Qualitative Results!

Inizializziamo rapidamente due array numpy ciascuno di dimensione 11 × 3 × 1, dove la prima dimensione rappresenta il numero di volumi che userai per testare il modello. La seconda dimensione rappresenta MSE e PSNR tra: immagini RM 7T di output predette e 7T di ground truth, output predetto e immagini RM 3T di input, immagini RM 3T di input e 7T di ground truth; la terza dimensione rappresenta il numero di canali che darai in input al modello!

mse= np.zeros([11,3,1])
psnr= np.zeros([11,3,1])
i=0 #for iterating over the slices of all the 11 volumes

Nella prossima parte della cella, itererai su tutti e 11 i volumi uno alla volta. In ogni iterazione, predirai ogni volume usando tutti e tre gli autoencoder e infine farai la media delle predizioni.

In ciascun volume, itererai sui canali del volume e calcolerai PSNR e MSE per i tre casi discussi sopra.

Poi, usando la libreria nibabel, salverai output predetto, input (3T) e ground truth (7T) come file in formato .nii: ciascuno degli 11 volumi composto da 207 slice.

Infine, salverai la matrice PSNR in un file di testo usando la libreria numpy.

Come indicato nel paper, fare la media degli output predetti aiuta a ridurre gli effetti del rumore preservando le caratteristiche locali nelle immagini ricostruite, migliorando così il PSNR rispetto agli output dei singoli decoder.

for j in range(11):
    decoded_imgs_1 = autoencoder_1.predict(ToPredict_images[i:i+207,:,:,:])
    decoded_imgs_2 = autoencoder_2.predict(ToPredict_images[i:i+207,:,:,:])
    decoded_imgs_3 = autoencoder_3.predict(ToPredict_images[i:i+207,:,:,:])
    decoded_imgs = np.mean( np.array([ decoded_imgs_1, decoded_imgs_2,decoded_imgs_3 ]), axis=0 )
    for channel in range(1):
        mse[j,0,channel]=  np.mean((ground_images[i:i+207,:,:,channel] - decoded_imgs[:,:,:,channel]) ** 2)
        psnr[j,0,channel] = 20 * math.log10( 1.0 / math.sqrt(mse[j,0,channel]))
        mse[j,1,channel]=  np.mean((ground_images[i:i+207,:,:,channel] - ToPredict_images[i:i+207,:,:,channel])** 2)
        psnr[j,1,channel] = 20 * math.log10( 1.0 / math.sqrt(mse[j,1,channel]))
        mse[j,2,channel] =  np.mean((ToPredict_images[i:i+207,:,:,channel] - decoded_imgs[:,:,:,channel]) ** 2)
        checklt = math.sqrt(mse[j,2,channel])
        psnr[j,2,channel] = 20 * math.log10( 1.0 / math.sqrt(mse[j,2,channel]))
    obj = nib.Nifti1Image(decoded_imgs, affine)
    string =str(j)+'_crossval1.nii'
    nib.save(obj, save + string)
    obj = nib.Nifti1Image(ground_images[i:i+207,:,:,:], affine)
    string =str(j)+'_ground_images_crossval1.nii'
    nib.save(obj, save + string)
    obj = nib.Nifti1Image(ToPredict_images[i:i+207,:,:,:], affine)
    string =str(j)+'_ToPredict_images_crossval1.nii'
    nib.save(obj, save + string)
    i=i+207


np.savetxt('psnr_all_slices.txt',psnr[:,:,0])

Se vuoi imparare di più su Python, segui il corso di DataCamp Introduction to Data Visualization with Matplotlib.

Dai un’occhiata al Keras Tutorial: Deep Learning in Python.

Argomenti
Python
Machine Learning
Deep Learning

Approfondisci Python e il Deep Learning

{count, plural, one {Corso} other {Corsi}}

Introduzione al Deep Learning con Keras

{hours} ore
46.4K
Inizia a creare modelli di deep learning con Keras.
Vedi i dettagliRight Arrow
Inizia Il Corso
Mostra altroRight Arrow
Correlato

{count, plural, one {blog} other {blogs}}

Tokenizzazione nel NLP: come funziona, sfide e casi d'uso

Guida al preprocessing NLP nel machine learning. Copriamo spaCy, i transformer di Hugging Face e come funziona la tokenizzazione in casi d'uso reali.
Abid Ali Awan's photo

Abid Ali Awan

10 min

{count, plural, one {blog} other {blogs}}

I 15 migliori server MCP remoti che ogni AI builder dovrebbe conoscere nel 2026

Scopri i 15 migliori server MCP remoti che stanno trasformando lo sviluppo AI nel 2026. Scopri come migliorano automazione, ragionamento, sicurezza e velocità dei workflow.
Abid Ali Awan's photo

Abid Ali Awan

15 min

{count, plural, one {blog} other {blogs}}

Che cos'è Snowflake? Guida per principianti alla piattaforma dati cloud

Esplora le basi di Snowflake, la piattaforma dati cloud. Scopri la sua architettura, le sue funzionalità e come integrarla nelle tue pipeline di dati.
Tim Lu's photo

Tim Lu

12 min

Mostra AltroMostra Altro