Super-résolution d’image avec un framework multi-décodeurs

Dans la partie Training Script, vous allez travailler sur un problème de super-résolution d’image en utilisant une architecture de deep learning inédite. L’objectif est d’apprendre une mise en correspondance non linéaire d’une IRM cérébrale à bas champ 3 Tesla vers une IRM à haut champ 7 Tesla. Dans cette partie Testing Script, vous utiliserez les poids entraînés en Partie 1 pour faire des prédictions sur des données inédites. Vous verrez aussi comment enregistrer des images 2D sous forme d’un volume combiné.

Le tutoriel est divisé en deux parties : la première couvre l’entraînement, la seconde le test.

Partie 1 : Training Script
Partie 2 : Testing Script

Remarque : Si vous souhaitez lire l’article, vous le trouverez ici.

Training Script

En bref, vous aborderez les points suivants dans la Partie 1 du tutoriel :

Commencez par importer les modules nécessaires pour entraîner votre modèle de deep learning,
Puis un aperçu du jeu de données IRM 3T et 7T,
Ensuite, vous définirez les initialiseurs et chargerez les jeux 3T et 7T ; au chargement, vous redimensionnerez aussi les images à la volée,
Puis vous préprosesserez les données chargées : conversion des listes train et test en matrices numpy, conversion des types en float32, remise à l’échelle min–max, reshape des tenseurs et découpe en 80 % entraînement et 20 % validation,
Vous créerez ensuite l’architecture 1-Encoder-3-Decoder avec connexions de fusion et décodeurs multiples,
Puis vous définirez la fonction de perte, créerez trois modèles distincts et les compilerez,
Enfin, vous entraînerez le modèle à connexions de fusion et multi-décodeurs, le testerez sur les données de validation et calculerez les résultats quantitatifs.

Dépendances des modules Python

Avant de suivre ce tutoriel, assurez-vous d’utiliser exactement les versions ci-dessous :

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

Remarque : Le modèle est entraîné sur une machine équipée d’un GPU Nvidia 1080 Ti, processeur Xeon e5 GeForce et 32 Go de RAM. Si vous utilisez Jupyter Notebook, ajoutez trois lignes pour définir l’ordre des périphériques CUDA et les GPU visibles via le module os.

Dans le code ci-dessous, vous définissez des variables d’environnement via os.environ dans le notebook. Il est recommandé de le faire avant d’initialiser Keras pour limiter TensorFlow (backend Keras) au premier GPU. Si la machine possède un GPU à l’index 0, utilisez 0 au lieu de 1. Vous pouvez le vérifier en lançant la commande nvidia-smi dans votre terminal.

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

Importation des modules

Commencez par importer les modules requis comme tensorflow, numpy et surtout keras, ainsi que les fonctions/couches nécessaires (Input, Conv2D, MaxPooling2D, etc.), puisque vous les utiliserez pour l’entraînement du modèle.
Pour lire des images au format NIfTI, importez également le module 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
Propulsé par

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)
Propulsé par

Comprendre le jeu de données IRM cérébrales 3T et 7T

Le jeu de données 3T et 7T se compose de volumes 3D ; chaque volume contient 207 coupes/images d’IRM cérébrales à différentes profondeurs. Chaque coupe mesure 173 x 173. Les images sont en niveaux de gris (monocanal). On compte 39 sujets au total, chaque sujet correspondant à une acquisition IRM d’un patient. Le format d’image n’est pas jpeg ou png, mais NIfTI. Vous verrez plus loin comment lire ce format.

Le jeu se compose d’images IRM en modalité T1, traditionnellement adaptées à l’évaluation des structures anatomiques. Dans ce tutoriel, vous travaillerez sur des IRM cérébrales 3T et 7T.

Le jeu de données est public et téléchargeable à cette adresse.

28 sujets sont utilisés pour l’entraînement et les 11 restants pour le test.

Définir les initialiseurs

Commençons par définir les dimensions des données. Vous redimensionnerez les images de 173x173 à 176x176 lors de la lecture. Ici, vous définissez aussi le répertoire des données, la taille de lot pour l’entraînement, le nombre de canaux, une couche Input(), des listes pour les matrices d’entraînement et de test, et enfin, pour la remise à l’échelle, vous chargez un fichier texte contenant les valeurs minimale et maximale du jeu IRM.

Remarque : Pour la remise à l’échelle, vous pouvez également prendre le min et le max de votre propre jeu de données si vous n’avez pas le fichier 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')Propulsé par

Chargement des données

Chargez ensuite les données IRM via la bibliothèque nibabel et redimensionnez les images de 173 x 173 à 176 x 176 en complétant par des zéros selon x et y.

Notez qu’au chargement d’un volume NIfTI, Nibabel ne lit pas immédiatement le tableau d’images ; il attend qu’on le lui demande, via la méthode get_data().

Comme vous souhaitez des coupes 2D (et non 3D), vous utiliserez les listes train et test initialisées plus haut ; à chaque volume lu, itérez sur les 207 coupes du volume 3D et ajoutez-les une par une à la liste.

folder = os.listdir(inp)
Propulsé par

Chargeons d’abord les données 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,:])
Propulsé par

Puis chargez les données 7T. Vous pouvez réutiliser la variable folder puisque le nombre de volumes 3T et 7T est identique.

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,:])
Propulsé par

Prétraitement des données

Les matrices train et test sont des listes ; convertissez-les en tableaux numpy.

Ensuite, convertissez leur type en float32 et remettez à l’échelle à l’aide de la stratégie min–max, pour l’input et le 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)
Propulsé par

Vérifions rapidement la forme de train_matrix (3T) et test_matrix (7T). Elles doivent contenir 28 x 207 = 5 796 images au total, chacune en 176 x 176.

train_matrix.shape
Propulsé par

(5796, 176, 176)
Propulsé par

test_matrix.shape
Propulsé par

(5796, 176, 176)
Propulsé par

Créez ensuite deux nouvelles variables : augmented_images (3T/input) et Haugmented_images (7T/ground truth) ayant la forme des matrices train et test. Ce seront des matrices 4D : première dimension = nombre total d’images, deuxième et troisième = dimensions de chaque image, dernière = nombre de canaux (ici 1).

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)])
Propulsé par

Passez ensuite sur toutes les images : à chaque itération, faites un reshape en 176 x 176 et affectez aux tenseurs augmented_images (3T/input) et 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)
Propulsé par

Il est ensuite crucial de partitionner les données. Pour une bonne généralisation, découpez en train/validation : 80 % pour l’entraînement et 20 % pour la validation.

Cela réduit aussi le risque de surapprentissage en validant sur des données non vues pendant l’entraînement.

Utilisez train_test_split de scikit-learn défini au début pour un découpage correct :

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)
Propulsé par

X_test = np.array(X_test)
y_test = np.array(y_test)
X_test = X_test.astype('float32')
y_test = y_test.astype('float32')
Propulsé par

X_train = np.array(X_train)
y_train = np.array(y_train)
X_train = X_train.astype('float32')
y_train = y_train.astype('float32')
Propulsé par

Le modèle : 1 encodeur – 3 décodeurs

Figure : rétropropagation d’autoencodeur sélectif

Image tirée de cet article.

Connexions de fusion

Définissez ensuite l’architecture proposée avec des blocs de couches de filtres successives suivies d’un max pooling dans la partie encodeur, comme ci-dessous. Pour reconstruire la taille d’origine en sortie, une couche d’upsampling est ajoutée dans chaque bloc de décodeur. Lors de l’upsampling, des artefacts peuvent apparaître par manque de détails dans l’entrée sous-échantillonnée du décodeur. Vous concaténez donc l’entrée du décodeur avec sa version mise à l’échelle issue de l’encodeur, pour fournir des détails pertinents et améliorer la reconstruction. L’ajout de connexions de fusion améliore sensiblement le PSNR (de l’ordre de 5 dB). Ce choix architectural s’inspire de cet article.

Multi-décodeurs

L’approche utilise un encodeur unique et plusieurs décodeurs avec une entrée monocanal. Trois couches convolutionnelles sont utilisées dans chaque bloc de l’encodeur et des trois décodeurs, suivies d’une couche de normalisation de lot pour la stabilité numérique.

Dans l’encodeur, le premier bloc convolutionnel comporte 32 filtres, puis le nombre de filtres est doublé à chaque bloc. Dans les décodeurs, le premier bloc a 256 filtres, puis ce nombre est divisé par deux à chaque bloc. La taille des filtres est 3 x 3 partout.

La fonction d’activation ReLU est utilisée dans toutes les couches sauf la dernière. Les données étant normalisées entre 0 et 1, une activation Sigmoid est utilisée en sortie finale.

On sait que les détails locaux à différentes échelles jouent un rôle clé en reconstruction d’image. L’architecture considère des images à plusieurs échelles via des couches hiérarchiques de sous-échantillonnage (max pooling) et suréchantillonnage (facteur 2) dans l’encodeur et les décodeurs. La représentation encodée après trois sous-échantillonnages projette l’entrée de haute dimension vers un espace latent.

Encodeur

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
Propulsé par

Décodeur 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
Propulsé par

Décodeur 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
Propulsé par

Décodeur 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
Propulsé par

Dans les 4 cellules ci-dessus, vous avez défini quatre fonctions : une pour l’encodeur et trois pour les décodeurs. Vous pourriez définir une fonction decoder() unique et l’appeler trois fois, mais pour la clarté et afin d’éviter des effets de hasard entre décodeurs, elles sont séparées.

Fonction de perte

Vous utiliserez une erreur quadratique moyenne en excluant les valeurs (pixels) égales à zéro dans le ground truth y_t et la prédiction y_p.

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)))))
Propulsé par

Définition et compilation du modèle

Commencez par appeler la fonction encoder en lui passant l’entrée. Comme l’architecture utilise des connexions de fusion, la fonction encoder renvoie les sorties de cinq couches convolution qui seront ensuite fusionnées séparément avec chaque décodeur.

conv5,conv4,conv3,conv2,conv1 = encoder(input_img)
Propulsé par

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())
Propulsé par

/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
Propulsé par

Entraîner le modèle

Vous enregistrerez les poids uniquement lorsque le peak signal-to-noise ratio sur la validation s’améliore. Pour cela, définissez une liste psnr_gray_channel avec une valeur initiale 1 (valeur factice largement dépassée dès les premières itérations).

Initialisez le taux d’apprentissage à 1e-3 et appliquez une décroissance : réduction de 10 % tous les 20 epochs.

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

Remarque : Le code qui suit doit idéalement être exécuté en un seul bloc, mais il est découpé ici pour faciliter la compréhension.

Le modèle est entraîné pendant 500 epochs. Le learning rate initial est 1e-3 comme défini. Vous enregistrerez aussi PSNR et MSE après chaque epoch entre 7T et la prédiction. Vous utiliserez K.set_value pour diminuer le taux d’apprentissage de 10 % toutes les 20 époques pour les trois modèles.

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)
Propulsé par

Brassez ensuite les paires 3T/7T pour éviter le surapprentissage lié à l’ordre des échantillons. Calculez le nombre de lots selon la batch_size définie, puis itérez sur 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):
Propulsé par

En plus des valeurs de PSNR, stockez aussi les pertes des trois autoencodeurs et de leurs décodeurs respectifs.

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')
Propulsé par

À chaque batch, vous souhaitez faire voir les 32 échantillons suivants au modèle : le code suivant s’en charge.

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)),:]
Propulsé par

Pour la stratégie d’autoencodeur minimum, après chaque epoch vous testez sur les données d’entraînement via test_on_batch de Keras, qui renvoie trois pertes. Affichez-les.

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))
Propulsé par

Six cas de figure sont possibles dans votre réseau :
loss_1 : perte de l’autoencodeur 1
loss_2 : perte de l’autoencodeur 2
loss_3 : perte de l’autoencodeur 3

loss_1 peut être plus petite que loss_2 et loss_3 : dans ce cas, vous n’entraînez que l’autoencodeur 1. Les parties encodeur des autoencodeurs 2 et 3 sont gelées (False) et seuls leurs décodeurs sont entraînés ; consignez toutes les pertes dans les fichiers.

loss_2 plus petite que loss_1 et loss_3 : vous n’entraînez que l’autoencodeur 2, et vous gelez les parties encodeur de 1 et 3, etc.

loss_3 plus petite que loss_1 et loss_2 : vous n’entraînez que l’autoencodeur 3, etc.

loss_1 égale à loss_2 : entraînez soit l’autoencodeur 1 soit le 2, gelez l’encodeur du 3 et celui non choisi entre 1 et 2, puis n’entraînez que leurs décodeurs ; enregistrez les pertes.

loss_2 égale à loss_3 : entraînez soit le 2 soit le 3, gelez l’encodeur du 1 et celui non choisi entre 2 et 3 ; enregistrez les pertes.

loss_3 égale à loss_1 : entraînez soit le 3 soit le 1, gelez l’encodeur du 2 et celui non choisi entre 3 et 1 ; enregistrez les pertes.

Figure : architecture du modèle

Image tirée de cet article.

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()
Propulsé par

Étape essentielle : comme vous avez mis certaines couches à False ci-dessus, il faut les remettre à True afin que tous les autoencodeurs soient utilisés par test_on_batch et qu’ils ne restent pas gelés.

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
Propulsé par

Vous enregistrerez les poids de trois manières, dont deux ci-dessous : d’une part après chaque 100 epochs, d’autre part après chaque epoch (écrasement).

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")
Propulsé par

Test sur les données de validation

Brassez d’abord les données de validation, puis appliquez les mêmes 6 conditions que précédemment. Quelle que soit la condition vraie, testez votre modèle avec l’autoencodeur correspondant.

Calculez ensuite MSE et PSNR entre decoded_imgs (prédictions) et ground truth, puis enregistrez-les dans des fichiers texte.

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))
Propulsé par

Ici, vous n’enregistrez les poids que si le PSNR entre la prédiction 7T et le ground truth (7T) dépasse la meilleure valeur précédente stockée dans 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)
Propulsé par

Enregistrer input, ground truth et prédictions : résultats quantitatifs

Définissez une matrice numpy temp de taille 176 x 528 pour sauvegarder 3 images côte à côte (176 x 176 chacune). Enregistrez l’une des images validation 3T, 7T et prédite, et multipliez par 255 car les images ont été normalisées entre 0 et 1.

Enfin, utilisez scipy pour enregistrer l’image.

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
Propulsé par

Fermez les fichiers PSNR et MSE à la fin.

myfile_valid_psnr_7T.close()
myfile_valid_mse_7T.close()
Propulsé par

Enfin, réduisez le taux d’apprentissage de 10 % de sa valeur courante toutes les 20 époques.

if jj % 20 ==0:
        learning_rate = learning_rate - learning_rate * 0.10
Propulsé par

Testing Script

Dans la partie 2 de ce tutoriel, vous verrez :

L’import des modules requis pour entraîner votre modèle,
Un rappel sur le jeu de données IRM 3T et 7T,
La définition des initialiseurs et le chargement du jeu de test 3T et 7T, avec redimensionnement à la volée,
Le prétraitement des données : conversion des listes en matrices numpy, passage en float32, remise à l’échelle min–max, reshape et (si besoin) découpe train/validation,
La création de l’architecture 1-Encoder-3-Decoder : connexions de fusion et multi-décodeurs,
La définition de la perte, la création de trois modèles et le chargement des poids entraînés,
Enfin, la prédiction sur des données inédites et l’enregistrement des résultats quantitatifs et qualitatifs. Vous apprendrez aussi à sauvegarder des images 2D en un volume combiné avec nibabel.

Dépendances des modules Python

Avant de commencer, assurez-vous d’utiliser les mêmes versions que ci-dessous :

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

Remarque : Le modèle est entraîné sur une machine avec GPU Nvidia 1080 Ti, processeur Xeon e5 GeForce et 32 Go de RAM. Sous Jupyter Notebook, ajoutez trois lignes de configuration CUDA via le module os.

Dans le code ci-dessous, vous définissez des variables d’environnement via os.environ. Faites-le avant d’initialiser Keras pour contraindre TensorFlow à utiliser le premier GPU. Si votre machine a un GPU en position 0, utilisez 0 au lieu de 1. Vérifiez via nvidia-smi.

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

Importation des modules

Importez d’abord tensorflow, numpy et surtout keras ainsi que les couches nécessaires (Input, Conv2D, MaxPooling2D, etc.), car vous les utiliserez pour l’entraînement du modèle.
Pour lire des images au format NIfTI, importez aussi 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
Propulsé par

Comprendre le jeu de données IRM cérébrales 3T et 7T

Le jeu 3T/7T comprend des volumes 3D ; chaque volume contient 207 coupes d’IRM cérébrales. Chaque coupe mesure 173 x 173. Les images sont monocanal en niveaux de gris. On compte 39 sujets au total. Le format est NIfTI (et non jpeg/png). Vous verrez plus loin comment lire ce format.

Il s’agit d’images T1, traditionnellement adaptées à l’étude de l’anatomie. Aujourd’hui, vous travaillerez sur des IRM cérébrales 3T et 7T.

Le jeu est public et disponible au téléchargement à cette source.

Définir les initialiseurs

Définissons d’abord les dimensions. Vous redimensionnerez de 173x173 à 176x176 lors de la lecture. Nous définissons aussi le répertoire des données, la taille des lots, le nombre de canaux, une couche Input(), des listes train/test et, pour la remise à l’échelle, le fichier texte contenant min et max du jeu IRM.

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

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

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

Charger les volumes de test

Chargez les données IRM avec nibabel et redimensionnez de 173 x 173 à 176 x 176 par padding en x et y.

Rappel : au chargement d’un volume NIfTI, Nibabel ne lit pas directement le tableau d’images ; utilisez get_data().

Comme vous voulez des coupes 2D, utilisez les listes train et test initialisées ; pour chaque volume lu, itérez sur les 207 coupes et ajoutez-les une à une à la liste.

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,:])
Propulsé par

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,:])
Propulsé par

Prétraitement des données

Les matrices de test 3T et 7T étant des listes, convertissez-les en tableaux numpy.

Convertissez ensuite le type en float32 et remettez à l’échelle input et 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)
Propulsé par

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

Créez ensuite deux variables ToPredict_images (test 3T/input) et ground_images (test 7T/ground truth) ayant la forme des matrices train/test : matrices 4D (nombre d’images, hauteur, largeur, canaux = 1).

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)])
Propulsé par

Puis itérez sur toutes les images : à chaque fois, faites un reshape en 176 x 176 et affectez à ToPredict_images (test 3T/input) et ground_images (test 7T/ground truth).

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

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

Le modèle

Figure : architecture du modèle

Image tirée de cet article.

Encodeur

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
Propulsé par

Décodeur

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
Propulsé par

Fonction de perte

Vous utiliserez une erreur quadratique moyenne en excluant les pixels nuls dans y_t (ground truth) et y_p (prédiction).

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)))))
Propulsé par

Définition du modèle et chargement des poids dans les trois autoencodeurs

Appelez d’abord la fonction encoder en lui passant l’entrée. Comme l’architecture utilise des connexions de fusion, l’encodeur retourne la sortie de cinq couches convolution qui seront fusionnées avec les sorties des trois décodeurs.

conv5,conv4,conv3,conv2,conv1 = encoder(input_img)
Propulsé par

Créez maintenant trois modèles et chargez les poids entraînés correspondants.

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

autoencoder_1.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE1.h5")
Propulsé par

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

autoencoder_2.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE2.h5")
Propulsé par

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

autoencoder_3.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE3.h5")
Propulsé par

Prédiction sur volumes de test : résultats quantitatifs et qualitatifs

Initialisons rapidement deux tableaux numpy de taille 11 x 3 x 1 : la première dimension correspond au nombre de volumes testés, la seconde aux métriques MSE et PSNR pour : (1) sortie prédite 7T vs. ground truth 7T, (2) input 3T vs. ground truth 7T, (3) input 3T vs. sortie prédite ; la troisième dimension correspond au nombre de canaux en entrée du modèle.

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

Dans la suite, vous itérerez sur les 11 volumes. À chaque itération, vous prédirez chaque volume avec les trois autoencodeurs puis ferez la moyenne des prédictions.

Pour chaque volume, vous parcourrez les canaux et calculerez PSNR et MSE pour les trois cas vus plus haut.

Avec nibabel, vous enregistrerez la sortie prédite, l’entrée (3T) et le ground truth (7T) au format .nii : chacun des 11 volumes comportant 207 coupes.

Enfin, vous enregistrerez la matrice PSNR dans un fichier texte via numpy.

Comme indiqué dans l’article, la moyenne des sorties prédites réduit le bruit tout en préservant les caractéristiques locales dans les images reconstruites, ce qui améliore le PSNR par rapport aux sorties individuelles des décodeurs.

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])
Propulsé par

Pour aller plus loin en Python, suivez le cours de DataCamp : Introduction à la visualisation de données avec Matplotlib.

Découvrez aussi le Tutoriel Keras : Deep Learning en Python.

Sujets

Python Apprentissage automatique Apprentissage profond

Training Script

Dépendances des modules Python

Importation des modules

Comprendre le jeu de données IRM cérébrales 3T et 7T

Définir les initialiseurs

Chargement des données

Prétraitement des données

Le modèle : 1 encodeur – 3 décodeurs

Connexions de fusion

Multi-décodeurs

Encodeur

Décodeur 1

Décodeur 2

Décodeur 3

Fonction de perte

Définition et compilation du modèle

Entraîner le modèle

Test sur les données de validation

Enregistrer input, ground truth et prédictions : résultats quantitatifs

Testing Script

Dépendances des modules Python

Importation des modules

Comprendre le jeu de données IRM cérébrales 3T et 7T

Définir les initialiseurs

Charger les volumes de test

Prétraitement des données

Le modèle

Encodeur

Décodeur

Fonction de perte

Définition du modèle et chargement des poids dans les trois autoencodeurs

Prédiction sur volumes de test : résultats quantitatifs et qualitatifs

Sujets

Approfondir Python et le deep learning

Contenus associés

Approfondir Python et le deep learning