> 译自：https://medium.com/coinmonks/visualizing-brain-imaging-data-fmri-with-python-e1d0358d9dba
## 简介
大脑是人类目前所知的最复杂的器官，为了很好的了解大脑这个器官，我们做了很多努力，核磁共振成像(Magnetic Resonance Image,MRI)技术就是其中的重要突破，通过MRI的方式，我们可以获得大脑的一些数据。
近年来，随着机器学习的兴起，医学数据与机器学习结合使用的情况越来越多，而要有效的使用好医学数据，其前提就是处理好这些数据，本文内容会重点介绍如何使用Python来处理与分析这些脑成像数据，不会涉及过多医学知识。
## sMRI与fMRI
脑成像相关的数据可以去SPM网站中下载，SPM的含义是统计参数映射(Statistical Paramtric Mapping)，MRI生成的数据其实就是一种参数映射数据，当然，更加方便的是在工作公众号中回复data3获得相应的数据与jupyter代码文件。
下载后，其中有4个文件，其中README开头的为描述文件，fM00223为功能性核磁共振(funciton MRI,fMRI)成像数据,sM00223为结构性核磁共振(structural MRI, sMRI)成像数据。通过描述文件可知，这些数据是一个人躺在MRI机器上听「双音节词」时大脑的成像数据。
为了方便理解后面数据处理的内容，有必要理解sMRI与fMRI是什么以及两者的差异。
### 结构性核磁共振(sMRI)
因为人的体内存有大量的水分子，而水分子中还有氢原子，sMRI其实就是利用氢原子来成像，这意味着人身体中的内脏、软组织等含有高水分与脂肪的器官会被清楚的扫描出来，而大脑就是这样的一个器官，通过sMRI可以清晰的看到大脑中的密集结构与大量细节，但sMRI的成像无法观察到大脑的运动情况，即无法判断那些部位目前是比较活跃的，只能给出大脑的结构细节。
如下图，科学家利用sMRI对人体腹腔进行成像，从图可以看出，腹腔的结构很明显。
### 功能性核磁共振(fMRI)
为了弥补sMRI的缺陷，fMRI应运而出，fMRI主要通过血氧浓度水平依赖(Blood-oxygen-level dependent，BOLD)来成像，一个器官的某个部位活跃，此时这个器官的这个部位就需要消耗更多的氧气，以此为依据，来进行成像。
通过fMRI的方式，我们可以很好的判断出大脑此时那些区域是活跃的，但这种活动并不等同于神经活动，但fMRI也有缺陷，即它的成像会损失大量的细节。
如下图，科学家通过fMRI，利用BOLD探索大脑活跃区域。
但从图中可以看出，大脑的细节几乎看不清晰，所以目前常规的方式是使用sMRI对大脑结构进行成像，而fMRI对大脑活跃区域进行成像。
## sMRI数据可视化处理
通常，神经影像文件都会以相应的规律将数据存储在固定文件格式的文件中，我们可以通过NiBabel库来读/写常见的神经影像文件中的数据。
sMRI对应的数据在sM00223文件夹下，进入文件夹可以发现有两种不同文件格式的数据，分别是.img与.hdr，这其实是医学影像领域常见的格式，主要用于「分析」，其中，.img作为数据文件，包含二进制的图像资料，而.hdr作为头文件，包含图像的元数据，但这两种格式现在逐渐被.nifti格式代替，这是因为.hdr头文件难以完全真实反映元数据。
使用NiBabel将sMRI获得的数据载入：
import nibabel
\# sM0223对应的文件
data\_path = './fMRI\_data/sM00223/'
files = os.listdir(data_path)
\# 读取其中的数据
data_all = \[\]
for data_file in files:
if data_file\[-3:\] == 'hdr':
data = nibabel.load(data\_path + data\_file).get_data()
\# 打印数据维数
print(data.shape)
\# -------结果
(256, 256, 54, 1)
可以看出，sMRI产生的是4维数据，但第4维其实没有包含任何信息，其说明了sMRI每次扫描会产生54个数据切片，每个切片对应图像的大小为256x256个体素。
> 体素：类似像素，像素表示的是二维图像的最小单位，而体素则用于三维成像空间。
为了方面可视化每个切片的数据，可以简单处理一下数据：
import numpy as np
data = np.rot90(data.squeeze(), 1)
print(data.shape)
\# -------结果
(256, 256, 54)
上述代码中，先通过numpy.squeeze()删除了数组中的单维条目，此时无用的第四维会被删除，接着使用numpy.rot90()将数组逆时针旋转了90度。
简单处理后，直接使用Matplotlib对每10个切片中的一个切片进行数据的绘制。
import matplotlib.pyplot as plt
\# 创建 6 个子图，不清楚其中概念，可以看本公众号关于Matplotlib的文章
fig, ax = plt.subplots(1, 6, figsize=\[18, 3\])
n = 0
slice = 0
for _ in range(6):
ax\[n\].imshow(data\[:, :, slice\], 'gray')
ax\[n\].set_xticks(\[\])
ax\[n\].set_yticks(\[\])
ax\[n\].set_title('Slice number: {}'.format(slice), color='r')
n += 1
slice += 10
fig.subplots_adjust(wspace=0, hspace=0)
plt.show()
这就是通过sMRI数据绘制出的大脑结构图了，其中第0层切片是最低的一个(接近脖子位置)，而第50层切片是最高的一个(接近头顶)，在第20层，可以看具有眼睛的切片。
## fMRI数据可视化处理
阅读README.txt可知fM00223数据集中，每张图像的大小为64x64个体素，采集的片数为64以及采集的卷数为96。知道了这些信息，就可以对fM00223数据集中的数据进行重构。
打开fM00223文件夹，可以发现确实正好有96个Hdr文件，这意味着每个文件包含了一个卷的所有片。
\# fMRI数据的基本信息
x_size = 64
y_size = 64
n_slice = 64
n_volumes = 96
\# 获得所有文件
data\_path = './fMRI\_data/fM00223/'
files = os.listdir(data_path)
\# 读取数据并进行重塑
data_all = \[\]
for data_file in files:
if data_file\[-3:\] == 'hdr':
data = nibabel.load(data\_path + data\_file).get_data()
# 将所有数据添加到list中，从而多了一个维度：时间维度
data\_all.append(data.reshape(x\_size, y\_size, n\_slice))
接着就可以通过Matplotlib可视化展示数据了，因为组成这些数据的是体素，是三维的图像，我们无法直接看到所有的数据，所有进行切片操作，通常会横切大脑从而获得冠状面(coronal)、横切面(transversal)和矢状面(sagittal)这3个平面，理解这三个概念很重要，如下图：
+ 冠状面为左，右方向将人体纵切为前后两部分的断面
+ 横切面指横向水平切割的平面
+ 矢状面是指将躯体纵断为左右两部分的解剖平面
看一下下面的代码：
\# 创建 3x6 的子图，大小为 18x11
fig, ax = plt.subplots(3, 6, figsize=\[18, 11\])
\# 组织数据以冠状平面进行可视化，第四维度为时间维度，这里取第一个时间点
coronal = np.transpose(data_all, \[1, 3, 2, 0\])
coronal = np.rot90(coronal, 1)
\# 组织数据以横切平面进行可视化
transversal = np.transpose(data_all, \[2, 1, 3, 0\])
transversal = np.rot90(transversal, 2)
\# 组织数据以矢状平面进行可视化
sagittal = np.transpose(data_all, \[2, 3, 1, 0\])
sagittal = np.rot90(sagittal, 1)
\# 可视化不同平面
n = 10
\# 对每个切片的6个切面进行操作
for i in range(6):
ax\[0\]\[i\].imshow(coronal\[:, :, n, 0\], cmap='gray')
ax\[0\]\[i\].set_xticks(\[\])
ax\[0\]\[i\].set_yticks(\[\])
if i == 0:
ax\[0\]\[i\].set_ylabel('coronal', fontsize=25, color='r')
n += 10
n = 5
for i in range(6):
ax\[1\]\[i\].imshow(transversal\[:, :, n, 0\], cmap='gray')
ax\[1\]\[i\].set_xticks(\[\])
ax\[1\]\[i\].set_yticks(\[\])
if i == 0:
ax\[1\]\[i\].set_ylabel('transversal', fontsize=25, color='r')
n += 10
n = 5
for i in range(6):
ax\[2\]\[i\].imshow(sagittal\[:, :, n, 0\], cmap='gray')
ax\[2\]\[i\].set_xticks(\[\])
ax\[2\]\[i\].set_yticks(\[\])
if i == 0:
ax\[2\]\[i\].set_ylabel('sagittal', fontsize=25, color='r')
n += 10
fig.subplots_adjust(wspace=0, hspace=0)
plt.show()
## 尾
感谢阅读，如果喜欢，点好看关注，如果有很想跟我说的话，欢迎直接跟我对话，个人微信：sighblue。
