3D-Volumendaten mit Matplotlib betrachten

Die meisten kennen Bilddaten, die mit normalen Kameras aufgenommen wurden („naturelle Bilder“ in der Fachliteratur), aber auch Aufnahmen mit Spezialgeräten wie Mikroskopen oder Teleskopen. Wenn du in Python mit Bildern arbeitest, ist die gebräuchlichste Art der Darstellung die Funktion imshow aus Matplotlib, der beliebtesten Plotting-Bibliothek für Python.

In diesem Tutorial zeigen wir dir, wie du diese Funktion erweiterst, um 3D-Volumendaten darzustellen – also einen Stapel von Bildern, die zusammen eine 3D-Struktur beschreiben. Beispiele: Magnetresonanztomographie (MRT) und Computertomographie (CT) erfassen 3D-Strukturen im Körper; Röntgen-Mikrotomographie zeigt 3D-Strukturen in Materialien wie Glas oder Metalllegierungen; und Lichtscheibenmikroskope messen fluoreszierende Partikel in biologischem Gewebe.

Wir zeigen dir, wie du einen MRT-Datensatz herunterlädst und die Schnitte mit matplotlib visualisierst. Du lernst:

• Daten mit Hilfe der Bibliothek nibabel zu laden, die einen Reader für das NIfTI-Dateiformat bereitstellt. Wenn du diesen Schritt überspringen und direkt plotten willst, erfährst du außerdem, wie du die Daten über scikit-image lädst.
• Erste Schritte zu einem voll funktionsfähigen Slice-Viewer mit der Event-Handler-API von matplotlib, und
• wie du deinen Code umschreibst, um dein neues Tool zu nutzen.
• Zum Schluss räumen wir den Workspace auf, wenn alles erledigt ist.

Vorab die Basics: So zeigst du Bilder mit Matplotlibs imshow an.

Wenn du mit Jupyter Notebook arbeitest, aktiviere zuerst den interaktiven Matplotlib-Modus:

%matplotlib notebook

Jetzt kannst du matplotlib importieren und Daten darstellen. Wir laden Beispieldaten aus dem data-Modul der scikit-image-Bibliothek:

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

Hinweis: Im interaktiven Notebook-Modus bleibt die geöffnete Figure aktiv, bis du sie über das Powertsymbol oben rechts deaktivierst. Mach das jeweils, bevor du mit dem nächsten Plot weitermachst.

Diese Bilder sind zweidimensional (2D), denn sie liegen entlang zweier Dimensionen: x und y; in NumPy-Sprache: Zeilen und Spalten oder r und c.

Manche Bilder sind 3D, sie haben eine zusätzliche Tiefen-Dimension (z bzw. Ebenen). Dazu gehören MR-Aufnahmen und die serielle Section Transmission Electron Microscopy (ssTEM), bei der die Probe dünn geschnitten (wie Salami) und jeder Schnitt separat aufgenommen wird.

Um solche Bilder in matplotlib zu betrachten, wählen wir einen Schnitt und zeigen nur diesen an. Probieren wir das mit frei verfüglichen MRT-Daten aus dem Netz.

Zwischenspiel: Daten beschaffen

Wir laden einen Datensatz aus Buchel und Friston, Cortical Interactions Evaluated with Structural Equation Modelling and fMRI (1997). Zuerst erstellen wir ein temporäres Verzeichnis für den Download. Denk daran, es nach der Analyse zu löschen! Wenn du den Datensatz behalten willst, ändere d auf einen dauerhaften Speicherort deiner Wahl.

import tempfile

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

Jetzt laden wir die Daten herunter:

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

Und entpacken sie aus der zip-Datei in unser temporäres Verzeichnis:

import zipfile

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

Wenn du dir den Inhalt ansiehst, findest du eine Reihe von ‘.hdr’- und ‘.img’-Dateien.

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

Das ist das NIfTI-Dateiformat, dafür brauchen wir einen Reader. Zum Glück bietet die hervorragende Bibliothek nibabel genau das. Installiere sie mit conda install -c conda-forge nibabel oder pip install nibabel, und dann:

import nibabel

Jetzt können wir das Bild einlesen und mit der Methode .get_data() ein NumPy-Array für die Ansicht holen:

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

Tipp: Wenn du direkt mit dem Plotten der MRT-Daten weitermachen willst, führe Folgendes aus:

from skimage import io

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

Zurück zum Plotten

Schauen wir uns nun einen Schnitt im Array an:

plt.imshow(struct_arr[75])

Uff! Das wirkt ziemlich gestaucht. Grund: Die Auflösung entlang der vertikalen Achse ist in vielen MRTs anders als entlang der horizontalen Achsen. Das lässt sich über den Parameter aspect in imshow korrigieren:

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

Damit es einfacher wird, transponieren wir die Daten und betrachten nur horizontale Schnitte, die kein solches Feintuning brauchen.

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

Hübsch! Um einen anderen Schnitt oder eine andere Achse zu sehen, brauchen wir allerdings wieder einen Aufruf von imshow:

plt.imshow(struct_arr2[5])

Diese ständigen Aufrufe werden schnell mühsam. Lange habe ich 3D-Volumen mit Tools außerhalb von Python betrachtet, z. B. mit ITK-SNAP. Aber es ist erstaunlich einfach, dem Matplotlib-Viewer 3D-Scrollen beizubringen! So kannst du 3D-Daten direkt in Python erkunden, ohne ständig zwischen Exploration und Analyse den Kontext zu wechseln.

Der Schlüssel ist die matplotlib-Event-Handler-API. Damit definierst du Aktionen auf dem Plot — inklusive Ändern der dargestellten Daten! — als Reaktion auf bestimmte Tastendrücke oder Mausklicks.

In unserem Fall binden wir die Tasten J und K an „vorheriger Schnitt“ und „nächster Schnitt“:

def previous_slice():
    pass

def next_slice():
    pass

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

Einfach genug! Jetzt müssen wir diese Aktionen implementieren und der Figure mitteilen, dass sie Tastendrücke mit process_key verarbeiten soll. Letzteres ist simpel: Wir nutzen die Canvas-Methode mpl_connect:

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

Die komplette Doku zu mpl_connect findest du hier, inklusive weiterer Events (z. B. Mausklicks).

Nach kurzem Stöbern findet man heraus, dass imshow ein AxesImage-Objekt zurückgibt, das „im“ Matplotlib-Axes-Objekt lebt – genauer in dessen Attribut .images, wo gezeichnet wird. Dieses Objekt hat eine praktische Methode set_array, die die angezeigten Bilddaten austauscht. Wir müssen also nur:

• einen beliebigen Index plotten und diesen speichern, z. B. als Laufzeitattribut am Axes-Objekt.
• die Funktionen next_slice und previous_slice schreiben, die den Index ändern und mit set_array den entsprechenden Schnitt des 3D-Volumens setzen.
• mit der Canvas-Methode draw die Figure neu zeichnen.

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

Probieren wir es aus!

multi_slice_viewer(struct_arr2)

Läuft! Allerdings wirst du merken, dass das Scrollen nach oben mit K auch die horizontale Skala staucht. Hm? (Passiert nur, wenn sich die Maus über dem Bild befindet.)

Hintergrund: Zusätzliche Event-Handler werden in Matplotlib einfach übereinandergelegt. In diesem Fall ist K ein eingebautes Shortcut, um die x-Achse auf log-Skalierung zu stellen. Wenn wir K exklusiv nutzen wollen, müssen wir es aus den Standard-Keymaps von matplotlib entfernen. Diese liegen als Listen im Dictionary plt.rcParams, dem Speicher für systemweite Default-Einstellungen:

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

Das heißt: Das Drücken einer der Tasten in der Liste (also <key1> oder <key2>) führt <command> aus.

Daher schreiben wir eine Helferfunktion, die die von uns benötigten Tasten überall aus diesem Dictionary entfernt. (Diese Funktion gibt es in matplotlib noch nicht, wäre aber sicher ein willkommener Beitrag!)

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)

Ein voll funktionsfähiger Slice-Viewer

Ok, schreiben wir unsere Funktion um und nutzen das neue 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])

Jetzt sollten wir alle Schnitte unseres MRT-Volumens ohne Störungen durch die Default-Keymap betrachten können!

multi_slice_viewer(struct_arr2)

Ein schönes Extra: Diese Methode funktioniert mit jedem Matplotlib-Backend! Auch in der IPython-Konsole bekommst du die gleiche Interaktion wie im Browser. Und ebenso in einer Qt- oder Tkinter-App mit eingebettetem Matplotlib-Plot. Dieses einfache Tool ermöglicht dir also, immer komplexere Anwendungen rund um Matplotlibs Visualisierungsmöglichkeiten zu bauen.

Zum Schluss

Vergiss nicht, aufzuräumen und das temporäre Verzeichnis zu löschen (falls du eines erstellt hast):

import shutil

# Remove the temporary directory
shutil.rmtree(d)