visualizando dados volumétricos 3D com Matplotlib

Muita gente já está familiarizada com dados de imagem capturados por câmeras comuns (frequentemente chamados de “imagens naturais” na literatura científica), mas também por instrumentos especializados, como microscópios ou telescópios. Ao trabalhar com imagens em Python, a forma mais comum de exibi-las é usando a função imshow do Matplotlib, a biblioteca de gráficos mais popular do Python.

Neste tutorial, vamos mostrar como estender essa função para exibir dados volumétricos 3D — pense nisso como uma pilha de imagens que, juntas, descrevem uma estrutura 3D. Por exemplo, exames de ressonância magnética (MRI) e tomografia computadorizada (CT) medem a estrutura 3D dentro do corpo humano; microtomografia por raios X mede a estrutura 3D dentro de materiais como vidro ou ligas metálicas; e microscópios de lâmina de luz medem partículas fluorescentes em tecidos biológicos.

Vamos demonstrar como baixar um conjunto de dados de MRI e exibir os cortes usando matplotlib. Você vai aprender a:

• Obter seus dados com a ajuda da biblioteca nibabel, que oferece um leitor para arquivos no formato NIfTI. Para quem quiser pular essa etapa e ir direto ao plot, este tutorial também mostra como obter os dados com scikit-image.
• Dar os primeiros passos rumo a um visualizador de cortes totalmente funcional com a API de manipuladores de evento do matplotlib, e
• Reescrever parte do seu código para aproveitar a nova ferramenta.
• Por fim, ver como limpar seu ambiente de trabalho quando terminar.

Mas antes, vamos ao básico: como exibir imagens com o imshow do Matplotlib.

Você pode começar ativando o modo interativo do matplotlib se estiver usando o Jupyter Notebook:

%matplotlib notebook

Agora importe o matplotlib e exiba alguns dados. Vamos carregar alguns exemplos incluídos no módulo data da biblioteca 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');

Observação: ao executar o matplotlib no modo interativo do notebook, a figura aberta continua sendo a única figura ativa até você desativá-la usando o ícone de energia no canto superior direito. Faça isso antes de avançar para cada novo gráfico.

Essas imagens são chamadas de 2D (bidimensionais) porque estão dispostas em 2 dimensões: x e y, ou, na linguagem do NumPy, linhas e colunas (r e c).

Algumas imagens são 3D, pois têm uma dimensão adicional de profundidade (z, ou planos). Isso inclui imagens de ressonância magnética (MRI) e microscopia eletrônica de transmissão por seções seriadas (ssTEM), em que a amostra é fatiada finamente, como um salame, e cada fatia é imageada separadamente.

Para visualizar esse tipo de imagem no matplotlib, precisamos escolher um corte e exibir apenas esse corte. Vamos testar isso com um conjunto de dados de MRI disponível gratuitamente online.

Intervalo: baixando os dados

Vamos baixar um conjunto de dados descrito em Buchel e Friston, Cortical Interactions Evaluated with Structural Equation Modelling and fMRI (1997). Primeiro, criamos um diretório temporário para fazer o download. Lembre-se de apagá-lo quando terminar a análise! Se quiser manter esse conjunto de dados para uso futuro, mude d para um caminho permanente de sua preferência.

import tempfile

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

Agora, vamos baixar os dados:

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 extrair o conteúdo do arquivo zip para nosso diretório temporário:

import zipfile

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

Se você conferir o conteúdo do arquivo, vai encontrar vários arquivos ‘.hdr’ e ‘.img’.

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

Eles estão no formato NIfTI, e vamos precisar de um leitor para isso. Felizmente, a excelente biblioteca nibabel fornece esse leitor. Instale com conda install -c conda-forge nibabel ou pip install nibabel e então:

import nibabel

Agora podemos finalmente ler nossa imagem e usar o método .get_data() para obter um array NumPy para visualização:

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

Dica: se quiser ir direto para o plot do MRI, execute as linhas abaixo:

from skimage import io

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

Voltando aos gráficos

Vamos olhar um corte desse array:

plt.imshow(struct_arr[75])

Uou! Parece meio amassado! Isso acontece porque a resolução no eixo vertical em muitos MRIs não é a mesma dos eixos horizontais. Dá para corrigir passando o parâmetro aspect para a função imshow:

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

Mas, para facilitar, vamos simplesmente transpor os dados e olhar apenas para os cortes horizontais, que não precisam desse ajuste.

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

Bonito! Claro que, para ver outro corte ou um corte em outro eixo, precisamos fazer outra chamada a imshow:

plt.imshow(struct_arr2[5])

Ficar repetindo essas chamadas cansa rápido. Por muito tempo, eu visualizava volumes 3D usando ferramentas fora do Python, como o ITK-SNAP. Mas dá para adicionar facilmente recursos de “rolagem” 3D ao visualizador do matplotlib! Assim, conseguimos explorar dados 3D dentro do Python, minimizando a troca de contexto entre exploração e análise.

O segredo é usar a API de manipuladores de evento do matplotlib, que permite definir ações no gráfico — incluindo trocar os dados exibidos! — em resposta a teclas específicas ou cliques do mouse.

No nosso caso, vamos associar as teclas J e K do teclado a “corte anterior” e “próximo corte”:

def previous_slice():
    pass

def next_slice():
    pass

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

Simples! Claro que precisamos decidir como implementar essas ações e informar à figura que ela deve usar a função process_key para processar as teclas! A segunda parte é fácil: basta usar o método mpl_connect do canvas da figura:

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

A documentação completa do mpl_connect está aqui, incluindo outros tipos de eventos que você pode associar (como cliques do mouse).

Descobri com um pouco de exploração que imshow retorna um objeto AxesImage, que “mora” dentro do objeto Axes do matplotlib onde acontece toda a renderização, no atributo .images. E esse objeto oferece um prático método set_array para trocar os dados da imagem exibida! Então, tudo o que precisamos fazer é:

• plotar um índice arbitrário e armazená-lo, talvez como um atributo em tempo de execução no objeto Axes.
• fornecer funções next_slice e previous_slice que alterem o índice e usem set_array para definir o respectivo corte do volume 3D.
• usar o método draw do canvas da figura para redesenhar a figura com os novos dados.

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

Vamos testar!

multi_slice_viewer(struct_arr2)

Funcionou! Muito bom! Mas, se você testar aí, vai notar que rolar para cima com K também “amassa” a escala horizontal do gráfico. Ué? (Isso só acontece se o mouse estiver sobre a imagem.)

O que rola é que adicionar manipuladores de evento no Matplotlib simplesmente os empilha. Neste caso, K é um atalho nativo para mudar a escala do eixo x para logarítmica. Se quisermos usar K exclusivamente, precisamos removê-la dos mapas de teclas padrão do matplotlib. Eles vivem como listas no dicionário plt.rcParams, o repositório do matplotlib para configurações padrão do sistema:

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

onde pressionar qualquer tecla da lista (isto é, <key1> ou <key2>) faz com que <command> seja executado.

Assim, vamos escrever uma função auxiliar para remover as teclas que quisermos usar, onde quer que apareçam nesse dicionário. (Essa função ainda não existe no matplotlib, mas provavelmente seria uma contribuição bem-vinda!)

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)

um visualizador de cortes completo

Beleza, vamos reescrever nossa função para usar essa nova ferramenta:

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

Agora, devemos conseguir ver todos os cortes do nosso volume de MRI sem interferência chata do keymap padrão!

multi_slice_viewer(struct_arr2)

Uma vantagem desse método é que ele funciona em qualquer backend do matplotlib! Então, se você testar no console do IPython, ainda terá a mesma interação que no navegador! E o mesmo vale para um app Qt ou Tkinter incorporando um gráfico do matplotlib. Essa ferramenta simples permite construir aplicações cada vez mais robustas em cima dos recursos de visualização do matplotlib.

antes de ir

Não esqueça de fazer a limpeza e apagar o diretório temporário (se você criou um):

import shutil

# Remove the temporary directory
shutil.rmtree(d)