Visualización de datos volumétricos 3D con Matplotlib

Seguro que conoces los datos de imagen tomados con cámaras corrientes (en la literatura científica se llaman «imágenes naturales»), pero también con instrumentos especializados, como microscopios o telescopios. Al trabajar con imágenes en Python, la forma más habitual de mostrarlas es con la función imshow de Matplotlib, la librería de gráficos más popular de Python.

En este tutorial, vamos a ampliar esa función para visualizar datos volumétricos 3D, que puedes imaginar como una pila de imágenes. Juntas, describen una estructura en tres dimensiones. Por ejemplo, la resonancia magnética (MRI) y la tomografía computarizada (CT) miden la estructura 3D del interior del cuerpo humano; la microtomografía de rayos X mide la estructura 3D dentro de materiales como el vidrio o las aleaciones metálicas; y los microscopios de lámina de luz miden partículas fluorescentes dentro de tejidos biológicos.

Te mostraremos cómo descargar un conjunto de datos de MRI y visualizar los cortes con matplotlib. Aprenderás a:

• Obtener tus datos con la ayuda de la librería nibabel, que incluye un lector para archivos en formato NIfTI. Si prefieres saltarte este paso y ponerte a graficar directamente, este tutorial también te enseña a cargar los datos con scikit-image.
• Dar los primeros pasos hacia un visor de cortes totalmente funcional con la API de control de eventos de matplotlib, y
• Reescribir parte de tu código para aprovechar tu nueva herramienta.
• Por último, verás cómo dejar el entorno limpio cuando hayas terminado.

Pero antes, repasemos lo básico: cómo mostrar imágenes con imshow de Matplotlib.

Si trabajas con Jupyter Notebook, puedes activar el modo interactivo de matplotlib:

%matplotlib notebook

Ahora ya puedes importar matplotlib y mostrar algunos datos. Vamos a cargar datos de ejemplo incluidos en el módulo data de la librería 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: cuando ejecutas matplotlib en modo interactivo dentro del notebook, la figura abierta se mantiene como la única activa hasta que la desactives con el icono de encendido en la esquina superior derecha. Asegúrate de hacerlo antes de pasar a cada nueva gráfica.

Estas imágenes son bidimensionales o 2D porque se disponen en dos dimensiones: x e y o, en la jerga de NumPy, filas y columnas (r y c).

Algunas imágenes son 3D porque tienen una dimensión adicional de profundidad (z, o planos). Es el caso de la resonancia magnética (MRI) y de la microscopía electrónica de transmisión por secciones seriadas (ssTEM), donde la muestra se corta en lonchas finas, como un salami, y se captura cada loncha por separado.

Para ver este tipo de imágenes en matplotlib, tenemos que elegir un corte y mostrar solo ese corte. Probémoslo con datos de MRI gratuitos disponibles en línea.

Interludio: cómo obtener los datos

Vamos a descargar un conjunto de datos descrito en Buchel y Friston, Cortical Interactions Evaluated with Structural Equation Modelling and fMRI (1997). Primero, creamos un directorio temporal donde descargarlos. ¡Recuerda borrarlo cuando termines el análisis! Si quieres conservar este dataset para más adelante, cambia d por una ruta permanente a tu elección.

import tempfile

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

Ahora, descarguemos los datos:

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

Y extráelos del archivo zip a nuestro directorio temporal:

import zipfile

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

Si miras el contenido del archivo, encontrarás un montón de ficheros «.hdr» e «.img».

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

Están en formato NIfTI y necesitaremos un lector para abrirlos. Por suerte, la excelente librería nibabel proporciona uno. Asegúrate de instalarla con conda install -c conda-forge nibabel o pip install nibabel, y luego:

import nibabel

Ahora sí, podemos leer la imagen y usar el método .get_data() para obtener un array de NumPy que podamos visualizar:

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

Truco: si quieres pasar directamente a graficar los datos de la MRI, ejecuta estas líneas:

from skimage import io

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

Volvemos a graficar

Veamos ahora un corte de ese array:

plt.imshow(struct_arr[75])

¡Vaya! ¡Se ve bastante aplastado! Esto ocurre porque la resolución en el eje vertical de muchas MRIs no coincide con la de los ejes horizontales. Podemos corregirlo pasando el parámetro aspect a imshow:

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

Pero, para simplificar, vamos a transponer los datos y a mirar solo los cortes horizontales, que no requieren este ajuste.

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

¡Bonito! Por supuesto, si queremos ver otro corte, o uno en otro eje, necesitamos otra llamada a imshow:

plt.imshow(struct_arr2[5])

Hacer tantas llamadas se vuelve tedioso enseguida. Durante mucho tiempo veía volúmenes 3D con herramientas externas a Python, como ITK-SNAP. Pero resulta que es bastante fácil añadir capacidad de «desplazamiento» 3D al visor de matplotlib. Así podemos explorar datos 3D dentro de Python, minimizando los cambios de contexto entre exploración y análisis.

La clave es usar la API de control de eventos de matplotlib, que nos permite definir acciones sobre la gráfica —¡incluido cambiar los datos representados!— en respuesta a determinadas pulsaciones de teclas o clics del ratón.

En nuestro caso, vamos a asociar las teclas J y K del teclado a «corte anterior» y «corte siguiente»:

def previous_slice():
    pass

def next_slice():
    pass

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

¡Así de simple! Claro que tenemos que implementar esas acciones y decirle a la figura que use la función process_key para procesar las pulsaciones de teclado. Esto último es sencillo: basta con usar el método mpl_connect del lienzo de la figura:

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

La documentación completa de mpl_connect está aquí, con otros tipos de eventos que puedes vincular (como clics del ratón).

Tras un poco de exploración, descubrí que imshow devuelve un objeto AxesImage, que vive «dentro» del objeto Axes de matplotlib donde ocurre todo el dibujo, en su atributo .images. Y este objeto ofrece un práctico método set_array para sustituir los datos de imagen mostrados. Así que solo necesitamos:

• graficar un índice arbitrario y guardarlo, por ejemplo, como un atributo en tiempo de ejecución del objeto Axes.
• proporcionar funciones next_slice y previous_slice que cambien el índice y usen set_array para establecer el corte correspondiente del volumen 3D.
• usar el método draw del lienzo de la figura para redibujarla con los nuevos datos.

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

¡Probémoslo!

multi_slice_viewer(struct_arr2)

¡Funciona! ¡Genial! Pero si lo pruebas, verás que al desplazarte hacia arriba con K también se aplasta la escala horizontal de la gráfica. ¿Eh? (Solo ocurre si el ratón está sobre la imagen).

Lo que pasa es que al añadir controladores de eventos en Matplotlib, simplemente se apilan unos sobre otros. En este caso, K es un atajo predefinido para cambiar el eje x a escala logarítmica. Si queremos usar K exclusivamente, tenemos que quitarla de los mapas de teclas por defecto de matplotlib. Estos viven como listas en el diccionario plt.rcParams, que es el repositorio de ajustes predeterminados de matplotlib:

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

donde pulsar cualquiera de las teclas de la lista (es decir, <key1> o <key2>) ejecutará <command>.

Por tanto, escribiremos una función auxiliar que elimine de este diccionario las teclas que queramos reservar para nuestro uso dondequiera que aparezcan. (Esta función aún no existe en matplotlib, ¡pero sería una contribución bienvenida!).

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 visor de cortes totalmente funcional

Vale, reescribamos nuestra función para aprovechar esta nueva herramienta:

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

Ahora deberíamos poder ver todos los cortes del volumen de MRI sin interferencias molestas del mapa de teclas por defecto.

multi_slice_viewer(struct_arr2)

Una ventaja de este método es que funciona con cualquier backend de matplotlib. Así que, si lo pruebas en la consola de IPython, tendrás la misma interacción que en el navegador. Y lo mismo ocurre con una app de Qt o Tkinter que integre una gráfica de matplotlib. Esta herramienta sencilla te permite construir aplicaciones cada vez más complejas sobre las capacidades de visualización de matplotlib.

Antes de irte

No olvides dejar todo recogido y borrar el directorio temporal (si lo creaste):

import shutil

# Remove the temporary directory
shutil.rmtree(d)