3D-volumetrische data bekijken met Matplotlib

De meesten van jullie kennen beelddata, gemaakt met gewone camera’s (in de wetenschappelijke literatuur vaak “natuurlijke afbeeldingen” genoemd), maar ook met gespecialiseerde instrumenten zoals microscopen of telescopen. Als je in Python met afbeeldingen werkt, is de meest gebruikelijke manier om ze te tonen de functie imshow van Matplotlib, Python’s populairste plottingbibliotheek.

In deze tutorial laten we je zien hoe je deze functie uitbreidt om 3D-volumetrische data weer te geven, wat je kunt zien als een stapel afbeeldingen. Samen beschrijven ze een 3D-structuur. Zo meten MRI- en CT-scans de 3D-structuur in het menselijk lichaam; röntgenmicrotomografie meet de 3D-structuur in materialen zoals glas of metaallegeringen; en light-sheet-microscopen meten fluorescerende deeltjes in biologisch weefsel.

We demonstreren hoe je een MRI-dataset downloadt en de slices weergeeft met matplotlib. Je leert hoe je

• Je data ophaalt met behulp van de bibliotheek nibabel, die een reader biedt voor bestanden in het NIfTI-bestandsformaat. Voor wie deze stap wil overslaan en direct wil plotten, laat deze tutorial ook zien hoe je de data met scikit-image kunt ophalen.
• Vervolgens zet je de eerste stappen richting een volledig functionele slice viewer met de event handler-API van matplotlib, en
• Herschrijf je wat code om je nieuwe tool te gebruiken.
• Tot slot zie je hoe je je werkruimte kunt opruimen als je klaar bent.

Maar eerst even de basis: hoe je afbeeldingen weergeeft met Matplotlib’s imshow.

Je kunt beginnen met het inschakelen van de interactieve matplotlib-modus als je met Jupyter Notebook werkt:

%matplotlib notebook

Nu kun je matplotlib importeren en wat data weergeven. Je laadt voorbeelddata in die is inbegrepen in de data-module van de scikit-image-bibliotheek:

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

Let op: Als je matplotlib in de interactieve notebookmodus draait, blijft de geopende figuur de enige actieve figuur totdat je deze uitschakelt met het aan/uit-symbool rechtsboven in de figuur. Zorg dat je dat doet voordat je doorgaat vanaf elke plot.

Deze afbeeldingen heten 2-dimensionale of 2D-afbeeldingen omdat ze zijn uitgezet langs 2 dimensies: x en y, of, in NumPy-termen, rijen en kolommen, of r en c.

Sommige afbeeldingen zijn 3D, omdat ze een extra diepte-dimensie hebben (z, of vlakken). Daaronder vallen MRI en seriële section transmission electron microscopy (ssTEM), waarbij een monster in dunne plakjes wordt gesneden, zoals een salami, en elk plakje apart wordt afgebeeld.

Om dergelijke afbeeldingen in matplotlib te bekijken, moeten we een slice kiezen en alleen die slice weergeven. Laten we het uitproberen met vrij beschikbare MRI-data online.

Intermezzo: de data ophalen

We gaan een dataset downloaden beschreven in Buchel en Friston, Cortical Interactions Evaluated with Structural Equation Modelling and fMRI (1997). Eerst maken we een tijdelijke map om de data in te downloaden. Vergeet niet deze te verwijderen als we klaar zijn met de analyse! Wil je deze dataset bewaren voor later gebruik, verander dan d naar een meer permanente maplocatie naar keuze.

import tempfile

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

Laten we nu de data downloaden:

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

En uitpakken uit het zip-bestand naar onze tijdelijke map:

import zipfile

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

Als je naar de daadwerkelijke inhoud van het bestand kijkt, zie je een hoop ‘.hdr’- en ‘.img’-bestanden.

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

Dit zijn bestanden in het NIfTI-formaat, en we hebben een reader nodig om ze te openen. Gelukkig biedt de uitstekende nibabel-bibliotheek zo’n reader. Installeer deze met conda install -c conda-forge nibabel of pip install nibabel, en dan:

import nibabel

Nu kunnen we eindelijk onze afbeelding inlezen en met de methode .get_data() een NumPy-array verkrijgen om te bekijken:

# 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: als je direct verder wilt met het plotten van de MRI-data, voer dan de volgende codelijnen uit:

from skimage import io

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

Terug naar het plotten

Laten we nu naar een slice in die array kijken:

plt.imshow(struct_arr[75])

Whoa! Dat ziet er behoorlijk geplet uit! Dat komt omdat de resolutie langs de verticale as in veel MRI’s niet hetzelfde is als langs de horizontale assen. Dat kunnen we oplossen door de parameter aspect door te geven aan de functie imshow:

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

Maar om het eenvoudiger te maken, zullen we de data gewoon transponeren en alleen naar de horizontale slices kijken, die zulke aanpassingen niet nodig hebben.

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

Mooi! Natuurlijk hebben we voor een andere slice, of een slice langs een andere as, weer een extra aanroep van imshow nodig:

plt.imshow(struct_arr2[5])

Al die aanroepen worden al snel nogal omslachtig. Lange tijd bekeek ik 3D-volumes met tools buiten Python, zoals ITK-SNAP. Maar het blijkt vrij eenvoudig om 3D-“scroll”-mogelijkheden toe te voegen aan de matplotlib-viewer! Zo kunnen we 3D-data binnen Python verkennen, waardoor je minder vaak hoeft te wisselen tussen data-exploratie en data-analyse.

De sleutel is de matplotlib event handler-API, waarmee we acties kunnen definiëren voor de plot — inclusief het wijzigen van de gegevens van de plot! — als reactie op bepaalde toetsaanslagen of muisklikken.

In ons geval koppelen we de J- en K-toets op het toetsenbord aan “vorige slice” en “volgende slice”:

def previous_slice():
    pass

def next_slice():
    pass

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

Zo simpel is het! Natuurlijk moeten we nog uitvinden hoe we deze acties daadwerkelijk implementeren én we moeten de figuur vertellen dat hij de functie process_key moet gebruiken om toetsaanslagen te verwerken! Dat laatste is eenvoudig: we hoeven alleen de figure canvas-methode mpl_connect te gebruiken:

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

De volledige documentatie voor mpl_connect vind je hier, inclusief welke andere soorten events je kunt koppelen (zoals muisklikken).

Na een beetje speurwerk ontdekte ik dat imshow een AxesImage-object retourneert, dat “binnenin” het matplotlib-Axes-object leeft waar al het tekenen plaatsvindt, in het attribuut .images. En dit object biedt een handige methode set_array die de afbeeldingsdata die wordt getoond, vervangt! Dus het enige wat we hoeven te doen is:

• een willekeurige index plotten en die index opslaan, bijvoorbeeld als een extra runtime-attribuut op het Axes-object.
• de functies next_slice en previous_slice schrijven die de index aanpassen en met set_array de bijbehorende slice van het 3D-volume instellen.
• de figure canvas-methode draw gebruiken om de figuur met de nieuwe data te hertekenen.

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

Laten we het proberen!

multi_slice_viewer(struct_arr2)

Dit werkt! Mooi! Maar als je dit thuis probeert, merk je dat omhoog scrollen met K ook de horizontale schaal van de plot samendrukt. Huh? (Dit gebeurt alleen als je muis boven de afbeelding hangt.)

Wat er gebeurt, is dat het toevoegen van event handlers aan Matplotlib ze simpelweg boven op elkaar stapelt. In dit geval is K een ingebouwde sneltoets om de x-as op een logaritmische schaal te zetten. Als we K exclusief willen gebruiken, moeten we hem verwijderen uit de standaard key maps van matplotlib. Deze staan als lijsten in de dictionary plt.rcParams, de opslagplaats van matplotlib voor standaard systeeminstellingen:

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

waarbij het indrukken van een van de toetsen in de lijst (dus <key1> of <key2>) ervoor zorgt dat <command> wordt uitgevoerd.

Daarom schrijven we een hulpfunctie om toetsen die we zelf willen gebruiken, overal waar ze in deze dictionary voorkomen, te verwijderen. (Deze functie bestaat nog niet in matplotlib, maar zou waarschijnlijk een welkome bijdrage zijn!)

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)

Een volledig functionele slice viewer

Oké, laten we onze functie herschrijven om gebruik te maken van deze nieuwe tool:

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

Nu zouden we alle slices in ons MRI-volume moeten kunnen bekijken zonder lastige interferentie van de standaard keymap!

multi_slice_viewer(struct_arr2)

Een fijn kenmerk van deze methode is dat hij werkt op elk matplotlib-backend! Dus als je dit uitprobeert in de IPython-terminalconsole, krijg je nog steeds dezelfde interactie als in de browser! En hetzelfde geldt voor een Qt- of Tkinter-app die een matplotlib-plot insluit. Deze eenvoudige tool stelt je dus in staat om steeds complexere applicaties te bouwen rond de visualisatiemogelijkheden van matplotlib.

Voor je gaat

Vergeet niet om achter jezelf op te ruimen en de tijdelijke map te verwijderen (als je er een hebt gemaakt):

import shutil

# Remove the temporary directory
shutil.rmtree(d)