Visualizzare dati volumetrici 3D con Matplotlib

La maggior parte di voi ha familiarità con i dati di immagini scattate con fotocamere comuni (spesso chiamate “immagini naturali” in letteratura scientifica), ma anche con strumenti specializzati come microscopi o telescopi. Quando lavori con le immagini in Python, il modo più comune per visualizzarle è usare la funzione imshow di Matplotlib, la libreria di plotting più popolare di Python.

In questo tutorial, ti mostreremo come estendere questa funzione per visualizzare dati volumetrici 3D, che puoi immaginare come un insieme di immagini impilate. Insieme, descrivono una struttura tridimensionale. Per esempio, la risonanza magnetica (MRI) e la tomografia computerizzata (CT) misurano la struttura 3D all’interno del corpo umano; la microtomografia a raggi X misura la struttura 3D all’interno di materiali come il vetro o le leghe metalliche; e i microscopi a foglio di luce misurano particelle fluorescenti all’interno di tessuti biologici.

Ti mostreremo come scaricare un dataset MRI e visualizzarne le slice usando matplotlib. Imparerai a

• Ottenere i dati con l’aiuto della libreria nibabel, che fornisce un lettore per file in formato NIfTI. Per chi volesse saltare questo passaggio e iniziare subito a plottare, questo tutorial mostrerà come ottenere i dati con scikit-image.
• Fare i primi passi verso un visualizzatore di slice completamente funzionale con l’API dei gestori di eventi di matplotlib, e
• Infine, vedere come ripulire l’ambiente di lavoro quando hai finito.

Ma prima, vediamo alcune basi: come visualizzare immagini con imshow di Matplotlib.

Puoi iniziare abilitando la modalità interattiva di matplotlib se stai lavorando con Jupyter Notebook:

%matplotlib notebook

Ora puoi importare matplotlib e visualizzare alcuni dati. Caricherai alcuni dati di esempio inclusi nel modulo data della libreria 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');

Nota: quando esegui matplotlib in modalità notebook interattiva, la figura aperta rimane l’unica figura attiva finché non la disabiliti usando il simbolo di accensione in alto a destra della figura. Assicurati di farlo prima di passare oltre ogni grafico.

Queste immagini sono chiamate 2-dimensionali o 2D perché sono disposte lungo 2 dimensioni: x e y, o, nel gergo di NumPy, righe e colonne o r e c.

Alcune immagini sono 3D, nel senso che hanno una dimensione di profondità aggiuntiva (z, o piani). Tra queste ci sono la risonanza magnetica (MRI) e la microscopia elettronica a trasmissione a sezioni seriali (ssTEM), in cui un campione viene tagliato in fettine sottili, come un salame, e ciascuna fetta viene acquisita separatamente.

Per visualizzare tali immagini in matplotlib, dobbiamo scegliere una slice e mostrare solo quella. Proviamolo su alcuni dati MRI liberamente disponibili online.

Interludio: ottenere i dati

Scaricheremo un dataset descritto in Buchel e Friston, Cortical Interactions Evaluated with Structural Equation Modelling and fMRI (1997). Per prima cosa, creiamo una directory temporanea in cui scaricare i dati. Ricordiamoci di eliminarla quando abbiamo finito l’analisi! Se vuoi conservare questo dataset per usi futuri, cambia d con una posizione di directory più permanente a tua scelta.

import tempfile

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

Ora, scarichiamo i dati:

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

E estraiamoli dal file zip nella nostra directory temporanea:

import zipfile

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

Se guardi il contenuto effettivo del file, troverai un mucchio di file ‘.hdr’ e ‘.img’.

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

Questi sono nel formato NIfTI, e ci serve un lettore per aprirli. Fortunatamente, l’ottima libreria nibabel fornisce un lettore del genere. Assicurati di installarla con conda install -c conda-forge nibabel o pip install nibabel, e poi:

import nibabel

Ora possiamo finalmente leggere la nostra immagine e usare il metodo .get_data() per ottenere un array NumPy da visualizzare:

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

Tip: se vuoi passare direttamente al plotting dei dati MRI, esegui le seguenti righe di codice:

from skimage import io

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

Torniamo al plotting

Ora diamo un’occhiata a una slice in quell’array:

plt.imshow(struct_arr[75])

Whoa! Sembra piuttosto schiacciata! Questo perché la risoluzione lungo l’asse verticale in molte MRI non è la stessa che lungo gli assi orizzontali. Possiamo risolverlo passando il parametro aspect alla funzione imshow:

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

Ma, per semplificare, trasporremo i dati e guarderemo solo le slice orizzontali, che non richiedono questi aggiustamenti.

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

Bello! Ovviamente, per vedere un’altra slice o una slice lungo un asse diverso, serve un’altra chiamata a imshow:

plt.imshow(struct_arr2[5])

Tutte queste chiamate diventano piuttosto noiose molto in fretta. Per molto tempo ho visualizzato volumi 3D usando strumenti esterni a Python, come ITK-SNAP. Ma, a quanto pare, è abbastanza facile aggiungere capacità di “scrolling” 3D al viewer di matplotlib! Questo ci permette di esplorare dati 3D direttamente in Python, riducendo al minimo il bisogno di cambiare contesto tra esplorazione dei dati e analisi.

La chiave è usare l’API dei gestori di eventi di matplotlib, che ci permette di definire azioni da eseguire sul grafico — incluso cambiare i dati del grafico! — in risposta a determinati tasti premuti o clic del mouse.

Nel nostro caso, associamo i tasti J e K della tastiera a “slice precedente” e “slice successiva”:

def previous_slice():
    pass

def next_slice():
    pass

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

Abbastanza semplice! Ovviamente, dobbiamo capire come implementare davvero queste azioni e dobbiamo dire alla figura che deve usare la funzione process_key per gestire le pressioni della tastiera! Quest’ultimo passaggio è semplice: basta usare il metodo mpl_connect del canvas della figura:

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

Puoi trovare la documentazione completa di mpl_connect qui, inclusi gli altri tipi di eventi a cui puoi associare azioni (come i clic del mouse).

Mi è bastata un po’ di esplorazione per scoprire che imshow restituisce un oggetto AxesImage, che vive “dentro” l’oggetto Axes di matplotlib, dove avviene tutto il disegno, nel suo attributo .images. E questo oggetto fornisce un comodo metodo set_array che sostituisce i dati immagine visualizzati! Quindi, tutto ciò che dobbiamo fare è:

• plottare un indice arbitrario e memorizzarlo, magari come attributo a runtime aggiuntivo sull’oggetto Axes.
• fornire funzioni next_slice e previous_slice che cambiano l’indice e usano set_array per impostare la slice corrispondente del volume 3D.
• usare il metodo draw del canvas della figura per ridisegnare la figura con i nuovi dati.

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

Proviamolo!

multi_slice_viewer(struct_arr2)

Funziona! Ottimo! Ma, se lo provi a casa, noterai che scorrere verso l’alto con K schiaccia anche la scala orizzontale del grafico. Eh? (Questo succede solo se il mouse è sopra l’immagine.)

Quello che accade è che aggiungere gestori di eventi a Matplotlib semplicemente li accumula uno sopra l’altro. In questo caso, K è una scorciatoia da tastiera integrata per cambiare l’asse x in scala logaritmica. Se vogliamo usare K esclusivamente, dobbiamo rimuoverla dalle mappature di tasti predefinite di matplotlib. Queste vivono come liste nel dizionario plt.rcParams, che è il repository di matplotlib per le impostazioni predefinite di sistema:

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

dove premere uno qualsiasi dei tasti nella lista (cioè <key1> o <key2>) farà eseguire <command>.

Quindi, dobbiamo scrivere una funzione di supporto per rimuovere i tasti che vogliamo usare ovunque possano comparire in questo dizionario. (Questa funzione non esiste ancora in matplotlib, ma probabilmente sarebbe un contributo gradito!)

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 visualizzatore di slice completamente funzionale

Ok, riscriviamo la nostra funzione per sfruttare questo nuovo strumento:

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

Ora dovremmo riuscire a visualizzare tutte le slice del nostro volume MRI senza fastidiose interferenze dalla keymap predefinita!

multi_slice_viewer(struct_arr2)

Una caratteristica interessante di questo metodo è che funziona su qualsiasi backend di matplotlib! Quindi, se lo provi nella console terminale di IPython, otterrai la stessa interazione che nel browser! E lo stesso vale per un’app Qt o Tkinter che incorpora un grafico matplotlib. Questo semplice strumento ti permette quindi di costruire applicazioni sempre più complesse attorno alle capacità di visualizzazione di matplotlib.

Prima di andare

Non dimentichiamo di ripulire dopo di noi ed eliminare la directory temporanea (se l’hai creata):

import shutil

# Remove the temporary directory
shutil.rmtree(d)