使用 Matplotlib 查看三维体数据

大多数人都熟悉由普通相机拍摄的图像数据（学术文献中常称为“自然图像”），也熟悉由显微镜或望远镜等专业仪器采集的图像。在 Python 中处理图像时，最常见的显示方式是使用 Matplotlib（Python 最流行的绘图库）的 imshow 函数。

在本教程中，我们将向您展示如何扩展该函数以显示三维体数据，您可以将其视为一叠图像。它们共同描述一个三维结构。例如，磁共振成像（MRI）和计算机断层扫描（CT）测量人体内部的三维结构；X 射线显微断层扫描测量玻璃或金属合金等材料内部的三维结构；光片显微镜测量生物组织内的荧光颗粒。

我们将演示如何下载一个 MRI 数据集，并使用 matplotlib 显示其切片。您将学习如何：

• 借助 nibabel 库获取数据，该库提供 NIfTI 文件格式的读取器。若您想跳过此步骤直接开始绘图，本教程也将展示如何使用 scikit-image 获取数据。
• 接着，使用 matplotlib 的事件处理 API，迈出打造完整切片查看器的第一步；以及
• 重构部分代码以使用您的新工具。
• 最后，了解完成后如何清理工作空间。

但首先，让我们回顾一些基础：如何使用 Matplotlib 的 imshow 显示图像。

如果您在 Jupyter Notebook 中工作，可以先启用 matplotlib 交互模式：

%matplotlib notebook

现在可以导入 matplotlib 并显示一些数据。您将加载 scikit-image 库 data 模块中包含的示例数据：

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

注意：当您在交互式笔记本模式下运行 matplotlib 时，打开的图形会一直保持为唯一的活动图形，直到您通过图形右上角的电源符号将其禁用。在继续下一幅图之前，请务必先执行该操作。

这些图像称为二维（2D）图像，因为它们沿两个维度展开：x 和 y，或用 NumPy 的惯用语，行和列（r 与 c）。

有些图像是三维（3D）的，即具有额外的深度维度（z，或称为平面）。这包括磁共振成像（MRI）和串行切片透射电子显微镜（ssTEM），其中样本被像萨拉米香肠一样薄切，每片分别成像。

要在 matplotlib 中查看此类图像，我们必须选择一个切片，仅显示该切片。让我们在网上一些可自由获取的 MRI 数据上试试。

插曲：获取数据

我们将下载 Buchel 和 Friston 于 1997 年发表的 Cortical Interactions Evaluated with Structural Equation Modelling and fMRI 一文中描述的数据集。首先，我们创建一个临时目录以下载数据。完成分析后务必删除它！如果您想保留该数据集以备后用，请将 d 更改为您选择的更持久的目录位置。

import tempfile

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

现在，来下载数据：

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

并将其从 zip 文件解压到我们的临时目录：

import zipfile

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

如果您查看文件的实际内容，会发现一堆“.hdr”和“.img”文件。

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

这些采用 NIfTI 文件格式，我们需要一个读取器。幸运的是，优秀的 nibabel 库提供了这样的读取器。请务必使用 conda install -c conda-forge nibabel 或 pip install nibabel 安装它，然后：

import nibabel

现在，我们终于可以读取图像，并使用 .get_data() 方法获取用于查看的 NumPy 数组：

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

提示：如果您想直接继续绘制 MRI 数据，请执行以下代码：

from skimage import io

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

回到绘图

现在让我们看看该数组中的一个切片：

plt.imshow(struct_arr[75])

哇！看起来挺“挤”的！这是因为许多 MRI 在垂直轴上的分辨率与水平方向不同。我们可以通过向 imshow 函数传递 aspect 参数来修复：

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

不过，为了更省事，我们将直接转置数据，只查看不需要这种调整的水平切片。

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

很漂亮！当然，如果要查看其他切片，或沿不同轴的切片，我们需要再次调用 imshow：

plt.imshow(struct_arr2[5])

这些调用很快就会变得相当繁琐。很长一段时间里，我会使用 Python 之外的工具查看三维体数据，例如 ITK-SNAP。但事实证明，为 matplotlib 查看器添加三维“滚动”功能非常容易！这使我们可以在 Python 内部探索三维数据，最大限度减少在数据探索与数据分析之间切换的需要。

关键在于使用 matplotlib 的事件处理 API，它允许我们为特定的按键或鼠标点击定义在绘图上执行的操作——包括更改绘图数据！

在我们的例子中，让我们将键盘上的 J 和 K 键绑定为“上一张切片”和“下一张切片”：

def previous_slice():
    pass

def next_slice():
    pass

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

相当简单！当然，我们需要想办法实际实现这些操作，并且要告诉图形应使用 process_key 函数来处理键盘按键！后者很简单：只需使用图形画布方法 mpl_connect：

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

您可以在此处找到 mpl_connect 的完整文档，包括可以绑定的其他事件类型（如鼠标点击）。

我稍加探索后发现，imshow 返回一个 AxesImage 对象，它“位于”进行所有绘制的 matplotlib Axes 对象内部的 .images 属性中。该对象提供了方便的 set_array 方法，可替换正在显示的图像数据！因此，我们只需：

• 绘制任意索引，并存储该索引，也许作为 Axes 对象的额外运行时属性。
• 提供 next_slice 和 previous_slice 函数，修改索引并使用 set_array 设置三维体对应的切片。
• 使用图形画布的 draw 方法以新数据重绘图形。

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

试试看！

multi_slice_viewer(struct_arr2)

能用！不错！不过，如果您在本地试试，会注意到向上滚动（按 K）也会把图的水平比例压缩。嗯？（只有当鼠标位于图像上方时才会发生。）

原因是向 Matplotlib 添加事件处理程序时，它们只是被层层叠加。在这个例子中，K 是一个内置快捷键，用于将 x 轴改为对数刻度。如果我们想专用 K，就必须将其从 matplotlib 的默认按键映射中移除。这些映射以列表形式存放在 plt.rcParams 字典中，它是 matplotlib 的系统级默认设置库：

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

按下列表中的任意按键（即 <key1> 或 <key2>）都会触发执行 <command>。

因此，我们需要编写一个辅助函数，从该字典中任何位置移除我们想要占用的按键。（该函数在 matplotlib 中尚不存在，但大概会是个受欢迎的贡献！）

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)

一个功能完整的切片查看器

好，让我们重写函数以使用这个新工具：

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

现在，我们应该可以在 MRI 体数据中查看所有切片，而不会再受到默认按键映射的干扰了！

multi_slice_viewer(struct_arr2)

这种方法的一个优点是它适用于任何 matplotlib 后端！因此，即使您在 IPython 终端控制台中尝试，也能获得与浏览器中相同的交互体验！对嵌入 matplotlib 图的 Qt 或 Tkinter 应用同样如此。这个简单工具因此使您能够围绕 matplotlib 的可视化能力构建更复杂的应用程序。

临走之前

别忘了善后，删除临时目录（如果您创建了的话）：

import shutil

# Remove the temporary directory
shutil.rmtree(d)