Tutoriel : super-résolution d’image avec un framework multi-décodeurs

Dans ce tutoriel, vous implémenterez un article de recherche en imagerie médicale avec Python dans Keras.

Actualisé 19 sept. 2026 · 15 min lire

Explorer avec l’IA

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

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

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)

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')

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)

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

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

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)

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

(5796, 176, 176)

test_matrix.shape

(5796, 176, 176)

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)])

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)

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)

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')

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

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

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

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

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)))))

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)

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

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

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)

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):

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')

À 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)),:]

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))

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()

É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

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")

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))

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)

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

Fermez les fichiers PSNR et MSE à la fin.

myfile_valid_psnr_7T.close()
myfile_valid_mse_7T.close()

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

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

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

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))

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')

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

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

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)

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

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)])

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)

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

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

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

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)))))

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)

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")

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")

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

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])

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