使用Matplotlib查看3D体积数据 _Matplotlib

本文概述

插曲：获取数据…
…返回绘图
功能齐全的切片查看器
你走之前…
总结

你们中的大多数人都熟悉用普通相机拍摄的图像数据(在科学文献中通常称为” 自然图像” ), 但也使用诸如显微镜或望远镜的专用仪器。使用Python处理图片时, 最常见的显示方式是使用Matplotlib的imshow函数, 该函数是Python最受欢迎的绘图库。
在本教程中, 我们将向你展示如何扩展此功能以显示3D体积数据, 你可以将其视为一堆图像。他们一起描述了3D结构。例如, 磁共振成像(MRI)和计算机断层扫描(CT)扫描可测量人体内部的3D结构。 X射线显微照相术可测量玻璃或金属合金等材料内部的3D结构；和光片显微镜测量生物组织内部的荧光颗粒。
我们将演示如何下载MRI数据集并使用matplotlib显示切片。你将学习如何

在nibabel库的帮助下获取数据, 该库为NIfTI文件格式的文件提供了阅读器。对于那些可能希望跳过此步骤并直接开始绘图的人, 本教程将向你展示如何使用scikit-image获取数据。
接下来, 你将迈出第一步, 使用matplotlib的事件处理程序API实现功能齐全的切片查看器, 以及
你将重写一些代码以使用新工具。
最后, 你将看到完成后如何清理工作区。

但是首先, 让我们看一些基本知识：如何使用Matplotlib的imshow显示图像。
如果你正在使用Jupyter Notebook, 则可以先启用交互式matplotlib模式：

%matplotlib notebook

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

import matplotlib.pyplot as pltfrom skimage import dataastronaut = data.astronaut()ihc = data.immunohistochemistry()hubble = data.hubble_deep_field()# Initialize the subplot panels side by sidefig, ax = plt.subplots(nrows=1, ncols=3)# Show an image in each subplotax[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图像, 因为它们沿2个维度布置：x和y, 或者以NumPy的说法是行和列或r和c。
一些图像是3D的, 因为它们具有附加的深度尺寸(z或平面)。其中包括磁共振成像(MRI)和连续切片透射电子显微镜(ssTEM), 其中像萨拉米香肠一样将样品切成薄片, 并对每个切片分别成像。
要在matplotlib中查看此类图像, 我们必须选择一个切片, 然后仅显示该切片。让我们在线上免费获取一些MRI数据进行尝试。
插曲：获取数据… 我们将下载Buchel和Friston在” 用结构方程模型和fMRI评估皮质相互作用” (1997)中描述的数据集。首先, 我们创建一个用于下载数据的临时目录。完成分析后, 我们必须记住要删除它！如果要保留此数据集供以后使用, 请将d更改为你选择的更永久的目录位置。

import tempfile# Create a temporary directoryd = tempfile.mkdtemp()

现在, 让我们下载数据：

import os# Return the tail of the pathos.path.basename('http://google.com/attention.zip')

from urllib.request import urlretrieve# Define URLurl = 'http://www.fil.ion.ucl.ac.uk/spm/download/data/attention/attention.zip'# Retrieve the datafn, info = urlretrieve(url, os.path.join(d, 'attention.zip'))

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

import zipfile# Extract the contents into the temporary directory we created earlierzipfile.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 metadatastruct_arr = struct.get_data()

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

from skimage import iostruct_arr = io.imread("https://s3.amazonaws.com/assets.srcmini.com/blog_assets/attention-mri.tif")

…返回绘图现在让我们看一下该数组中的一个切片：

plt.imshow(struct_arr[75])

文章图片
哇！看起来很松软！这是因为在许多MRI中, 沿纵轴的分辨率与沿横轴的分辨率不同。我们可以通过将Aspect参数传递给imshow函数来解决此问题：

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

文章图片
但是, 为了使事情变得简单, 我们将转置数据并仅查看水平切片, 而无需进行此类操作。

struct_arr2 = struct_arr.Tplt.imshow(struct_arr2[34])

文章图片
【使用Matplotlib查看3D体积数据】漂亮！当然, 要查看另一个切片或沿不同轴的切片, 我们需要另一个调用以显示：

plt.imshow(struct_arr2[5])

文章图片
所有这些电话很快就变得很乏味。长期以来, 我将使用Python之外的工具(例如ITK-SNAP)查看3D体积。但是, 事实证明, 向matplotlib查看器添加3D” 滚动” 功能非常容易！这使我们能够在Python中探索3D数据, 从而最大限度地减少了在数据探索和数据分析之间切换上下文的需求。
关键是使用matplotlib事件处理程序API, 该API可让我们定义要在绘图上执行的操作-包括更改绘图的数据！ —响应特定的按键或鼠标单击。
在我们的情况下, 让我们将键盘上的J和K键绑定到” 上一个切片” 和” 下一个切片” ：

def previous_slice():passdef next_slice():passdef 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对象, 该对象位于.pimages属性中发生所有绘制的matplotlib Axes对象” 内部” 。这个对象提供了一个方便的set_array方法, 可以交换出正在显示的图像数据！因此, 我们需要做的是：

绘制任意索引, 然后将该索引存储为Axes对象上的其他运行时属性。
提供函数next_slice和previous_slice, 这些函数可更改索引并使用set_array设置3D体的相应切片。
使用图形画布绘制方法使用新数据重新绘制图形。

def multi_slice_viewer(volume):fig, ax = plt.subplots()ax.volume = volumeax.index = volume.shape[0] // 2ax.imshow(volume[ax.index])fig.canvas.mpl_connect('key_press_event', process_key)def process_key(event):fig = event.canvas.figureax = 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.volumeax.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.volumeax.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_setfor key in remove_list:keys.remove(key)

功能齐全的切片查看器好的, 让我们重写我们的功能以使用此新工具：

def multi_slice_viewer(volume):remove_keymap_conflicts({'j', 'k'})fig, ax = plt.subplots()ax.volume = volumeax.index = volume.shape[0] // 2ax.imshow(volume[ax.index])fig.canvas.mpl_connect('key_press_event', process_key)def process_key(event):fig = event.canvas.figureax = 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.volumeax.index = (ax.index - 1) % volume.shape[0]# wrap around using %ax.images[0].set_array(volume[ax.index])def next_slice(ax):volume = ax.volumeax.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 directoryshutil.rmtree(d)

总结恭喜！这是一段漫长的旅程, 但是你已经完成了本Matplotlib教程的结尾！
请确保继续关注本系列的第二篇文章, 在该文章中, 你将了解有关缩放比例图, 显示每幅图的切片位置的十字线以及鼠标交互性的更多信息。
如果你想进一步了解Python的数据可视化, 请考虑参加srcmini的Python数据可视化入门课程。