Visualiser des données volumétriques 3D avec Matplotlib

La plupart d’entre vous connaissent les images issues d’appareils photo classiques (souvent appelées « images naturelles » dans la littérature scientifique), mais aussi celles provenant d’instruments spécialisés, comme les microscopes ou les télescopes. En Python, la manière la plus courante de les afficher consiste à utiliser la fonction imshow de Matplotlib, la bibliothèque de visualisation la plus populaire en Python.

Dans ce tutoriel, nous allons vous montrer comment étendre cette fonction pour afficher des données volumétriques 3D, que vous pouvez voir comme une pile d’images. Ensemble, elles décrivent une structure en trois dimensions. Par exemple, l’imagerie par résonance magnétique (IRM) et la tomodensitométrie (CT) mesurent la structure 3D à l’intérieur du corps humain ; la microtomographie aux rayons X mesure la structure 3D de matériaux comme le verre ou les alliages métalliques ; et les microscopes à nappe de lumière mesurent des particules fluorescentes dans des tissus biologiques.

Nous vous montrerons comment télécharger un jeu de données IRM et afficher les coupes avec matplotlib. Vous apprendrez à :

• Récupérer vos données grâce à la bibliothèque nibabel, qui fournit un lecteur pour les fichiers au format NIfTI. Pour celles et ceux qui souhaitent aller plus vite et tracer directement, ce tutoriel vous montrera aussi comment obtenir les données avec scikit-image.
• Faire vos premiers pas vers un visualiseur de coupes pleinement fonctionnel avec l’API de gestion d’événements de matplotlib, puis
• Réécrire une partie de votre code pour tirer parti de votre nouvel outil.
• Enfin, voir comment nettoyer votre espace de travail une fois terminé.

Mais d’abord, revoyons les bases : comment afficher des images avec imshow de Matplotlib.

Vous pouvez commencer par activer le mode interactif de matplotlib si vous travaillez dans Jupyter Notebook :

%matplotlib notebook

Vous pouvez maintenant importer matplotlib et afficher quelques données. Nous allons charger des exemples inclus dans le module data de la bibliothèque scikit-image :

import matplotlib.pyplot as plt
from skimage import data

astronaut = data.astronaut()
ihc = data.immunohistochemistry()
hubble = data.hubble_deep_field()

# Initialize the subplot panels side by side
fig, ax = plt.subplots(nrows=1, ncols=3)

# Show an image in each subplot
ax[0].imshow(astronaut)
ax[0].set_title('Natural image')
ax[1].imshow(ihc)
ax[1].set_title('Microscopy image')
ax[2].imshow(hubble)
ax[2].set_title('Telescope image');

Remarque : lorsque vous exécutez matplotlib en mode notebook interactif, la figure ouverte reste la seule figure active jusqu’à ce que vous la désactiviez via l’icône d’alimentation en haut à droite de la figure. Pensez à le faire avant de passer à chaque nouveau tracé.

Ces images sont dites bidimensionnelles ou 2D car elles s’étendent le long de 2 dimensions : x et y, ou, en termes NumPy, lignes et colonnes (r et c).

Certaines images sont 3D, car elles possèdent une dimension de profondeur supplémentaire (z, ou plans). C’est le cas notamment de l’imagerie par résonance magnétique (IRM) et de la microscopie électronique en transmission sériée (ssTEM), où un échantillon est finement tranché, comme un saucisson, et chaque tranche est imagée séparément.

Pour afficher ce type d’images avec matplotlib, il faut choisir une coupe et n’afficher que celle-ci. Essayons sur des données IRM libres d’accès en ligne.

Intermède : récupérer les données

Nous allons télécharger un jeu de données décrit dans Buchel et Friston, Cortical Interactions Evaluated with Structural Equation Modelling and fMRI (1997). Commençons par créer un répertoire temporaire pour y télécharger les données. N’oubliez pas de le supprimer une fois l’analyse terminée ! Si vous souhaitez conserver ce jeu de données, remplacez d par un emplacement plus pérenne de votre choix.

import tempfile

# Create a temporary directory
d = tempfile.mkdtemp()

Passons au téléchargement des données :

import os

# Return the tail of the path
os.path.basename('http://google.com/attention.zip')

from urllib.request import urlretrieve

# Define URL
url = 'http://www.fil.ion.ucl.ac.uk/spm/download/data/attention/attention.zip'

# Retrieve the data
fn, info = urlretrieve(url, os.path.join(d, 'attention.zip'))

Puis extrayons-les du fichier zip vers notre répertoire temporaire :

import zipfile

# Extract the contents into the temporary directory we created earlier
zipfile.ZipFile(fn).extractall(path=d)

Si vous inspectez le contenu, vous trouverez une série de fichiers « .hdr » et « .img ».

# List first 10 files
[f.filename for f in zipfile.ZipFile(fn).filelist[:10]]

Ils sont au format NIfTI, pour lequel il nous faut un lecteur. Heureusement, l’excellente bibliothèque nibabel en fournit un. Assurez-vous de l’installer avec conda install -c conda-forge nibabel ou pip install nibabel, puis :

import nibabel

Nous pouvons enfin lire notre image et utiliser la méthode .get_data() pour obtenir un tableau NumPy à afficher :

# Read the image 
struct = nibabel.load(os.path.join(d, 'attention/structural/nsM00587_0002.hdr'))

# Get a plain NumPy array, without all the metadata
struct_arr = struct.get_data()

Astuce : si vous souhaitez passer directement au tracé des données IRM, exécutez les lignes suivantes :

from skimage import io

struct_arr = io.imread("https://s3.amazonaws.com/assets.datacamp.com/blog_assets/attention-mri.tif")

Retour aux tracés

Regardons maintenant une coupe dans ce tableau :

plt.imshow(struct_arr[75])

Oups ! L’image semble toute « écrasée » ! C’est parce que la résolution verticale de nombreuses IRM diffère de celle des axes horizontaux. On peut corriger cela en passant le paramètre aspect à la fonction imshow :

plt.imshow(struct_arr[75], aspect=0.5)

Pour nous simplifier la vie, nous allons simplement transposer les données et ne regarder que les coupes horizontales, qui n’ont pas besoin de ce type d’ajustement.

struct_arr2 = struct_arr.T
plt.imshow(struct_arr2[34])

Sympa ! Évidemment, pour visualiser une autre coupe, ou une coupe selon un autre axe, il faut refaire un appel à imshow :

plt.imshow(struct_arr2[5])

À force, ces appels deviennent vite fastidieux. Longtemps, j’ai visualisé des volumes 3D avec des outils hors Python, comme ITK-SNAP. Mais il est en fait assez simple d’ajouter un « défilement » 3D au visualiseur matplotlib ! Cela permet d’explorer des données 3D directement dans Python, en limitant les changements de contexte entre exploration et analyse.

La clé, c’est d’utiliser l’API de gestion d’événements de matplotlib, qui nous permet de définir des actions à réaliser sur le graphique — y compris changer les données affichées ! — en réponse à certaines frappes clavier ou clics de souris.

Dans notre cas, lions les touches J et K du clavier à « coupe précédente » et « coupe suivante » :

def previous_slice():
    pass

def next_slice():
    pass

def process_key(event):
    if event.key == 'j':
        previous_slice()
    elif event.key == 'k':
        next_slice()

Simple ! Évidemment, il faut implémenter ces actions et indiquer à la figure d’utiliser la fonction process_key pour traiter les frappes clavier ! Ce dernier point est simple : il suffit d’utiliser la méthode mpl_connect du canevas de la figure :

fig, ax = plt.subplots()
ax.imshow(struct_arr[..., 43])
fig.canvas.mpl_connect('key_press_event', process_key)

La documentation complète de mpl_connect est disponible ici, avec la liste des autres événements possibles (comme les clics de souris).

Après quelques explorations, j’ai découvert que imshow renvoie un objet AxesImage, qui vit « dans » l’objet Axes de matplotlib où tout le dessin s’effectue, via son attribut .images. Et cet objet propose une méthode pratique, set_array, pour remplacer l’image affichée ! Il nous suffit donc de :

• tracer un indice quelconque et stocker cet indice, par exemple comme attribut d’exécution supplémentaire sur l’objet Axes ;
• fournir des fonctions next_slice et previous_slice qui modifient l’indice et utilisent set_array pour afficher la coupe correspondante du volume 3D ;
• utiliser la méthode draw du canevas pour redessiner la figure avec les nouvelles données.

def multi_slice_viewer(volume):
    fig, ax = plt.subplots()
    ax.volume = volume
    ax.index = volume.shape[0] // 2
    ax.imshow(volume[ax.index])
    fig.canvas.mpl_connect('key_press_event', process_key)

def process_key(event):
    fig = event.canvas.figure
    ax = fig.axes[0]
    if event.key == 'j':
        previous_slice(ax)
    elif event.key == 'k':
        next_slice(ax)
    fig.canvas.draw()

def previous_slice(ax):
    """Go to the previous slice."""
    volume = ax.volume
    ax.index = (ax.index - 1) % volume.shape[0]  # wrap around using %
    ax.images[0].set_array(volume[ax.index])

def next_slice(ax):
    """Go to the next slice."""
    volume = ax.volume
    ax.index = (ax.index + 1) % volume.shape[0]
    ax.images[0].set_array(volume[ax.index])

Essayons !

multi_slice_viewer(struct_arr2)

Ça marche ! Super ! Mais si vous essayez chez vous, vous noterez qu’en faisant défiler vers le haut avec K, l’échelle horizontale du tracé se compresse aussi. Hein ? (Cela se produit uniquement si votre souris est au-dessus de l’image.)

En fait, ajouter des gestionnaires d’événements dans Matplotlib les empile simplement les uns sur les autres. Ici, K est un raccourci clavier intégré pour passer l’axe des x en échelle logarithmique. Si nous voulons réserver K à notre usage, nous devons le retirer des mappages clavier par défaut de matplotlib. Ceux-ci sont enregistrés sous forme de listes dans le dictionnaire plt.rcParams, le dépôt des paramètres système par défaut de matplotlib :

plt.rcParams['keymap.<command>'] = ['<key1>', '<key2>']

où l’appui sur n’importe quelle touche de la liste (c.-à-d. <key1> ou <key2>) déclenche l’exécution de <command>.

Nous allons donc écrire une fonction utilitaire pour supprimer, partout où elles apparaissent, les touches que nous voulons réaffecter dans ce dictionnaire. (Cette fonction n’existe pas encore dans matplotlib, mais serait sans doute une contribution bienvenue !)

def remove_keymap_conflicts(new_keys_set):
    for prop in plt.rcParams:
        if prop.startswith('keymap.'):
            keys = plt.rcParams[prop]
            remove_list = set(keys) & new_keys_set
            for key in remove_list:
                keys.remove(key)

Un visualiseur de coupes pleinement fonctionnel

Réécrivons notre fonction pour utiliser ce nouvel outil :

def multi_slice_viewer(volume):
    remove_keymap_conflicts({'j', 'k'})
    fig, ax = plt.subplots()
    ax.volume = volume
    ax.index = volume.shape[0] // 2
    ax.imshow(volume[ax.index])
    fig.canvas.mpl_connect('key_press_event', process_key)

def process_key(event):
    fig = event.canvas.figure
    ax = fig.axes[0]
    if event.key == 'j':
        previous_slice(ax)
    elif event.key == 'k':
        next_slice(ax)
    fig.canvas.draw()

def previous_slice(ax):
    volume = ax.volume
    ax.index = (ax.index - 1) % volume.shape[0]  # wrap around using %
    ax.images[0].set_array(volume[ax.index])

def next_slice(ax):
    volume = ax.volume
    ax.index = (ax.index + 1) % volume.shape[0]
    ax.images[0].set_array(volume[ax.index])

Nous devrions maintenant pouvoir parcourir toutes les coupes de notre volume IRM sans interférences gênantes des raccourcis par défaut !

multi_slice_viewer(struct_arr2)

Atout appréciable : cette méthode fonctionne avec n’importe quel backend matplotlib ! Si vous l’essayez dans la console IPython, vous obtiendrez la même interaction que dans le navigateur ! Et c’est aussi vrai pour une application Qt ou Tkinter qui intègre un tracé matplotlib. Cet outil simple vous permet donc de bâtir des applications toujours plus riches autour des capacités de visualisation de matplotlib.

Avant de partir

N’oublions pas de faire le ménage et de supprimer le répertoire temporaire (si vous en avez créé un) :

import shutil

# Remove the temporary directory
shutil.rmtree(d)