1. fMRI相关知识

1.1 DICOM简介

​ DICOM （Digital Imaging and Communications in Medicine） 即医学数字成像和通信，是医学图像和相关信息的国际标准（ISO 12052）。 DICOM 被广泛应用于放射医疗，心血管成像以及放射诊疗诊断设备（X射线，CT，核磁共振，超声等），并且在眼科和牙科等其它医学领域得到越来越深入广泛的应用。目前采用的标准是 DICOM3.0 ，每一张图像中都携带着大量的信息。

1.2 NlfTI

​ 医学成像数据的文件格式，如 DICOM 和 NIfTI ，虽然类似但具有不同的功能。 NIfTI 主要用于医学成像和神经科学研究领域，而 DICOM 是存储和分发医学图片的标准文件格式。NIfTI 格式是为了改进以前用于存储医学成像数据的标准 Analyze 格式的缺点而创建的。 除了实际的图片数据，NIfTI 文件还包括元数据，如图像的测量、像素大小和方向，以及可能在原始数据上执行的任何预处理的结果。

​ DICOM 文件通常以一系列文件的形式存在，每个文件代表一层扫描。这意味着一个完整的3D扫描可能由数百个单独的文件组成。将它们转换为单一的 NIfTI 文件可以简化文件管理，使得加载和处理数据更加高效。有许多只能处理 NIfTI 文件的医学成像软件程序和分析工具，因此转换 DICOM 文件是一项常规任务。使用.nii.gz格式可以有效减小文件大小，便于存储和传输，同时保持数据的完整性和可用性。可以使用多种工具和软件包将 DICOM 转换为 NIfTI 。

  • dcm2niix - An improved variant of dcm2nii, dcm2niix can process more intricate DICOM data and produce JSON sidecar files.

    安装： conda install -c conda-forge dcm2niix

    常用命令： dcm2niix -z y -f %i_%d_%s -o outdir dcmdir

    1
    2
    3
    4
    5
    6
    7
    8
    -f 输出文件名 (default: -f %f%p%t%s)
    %i：插入病人ID（DICOM标签0010,0020）。例如，“myName％i”的输出文件名将把病人ID命名为“ID123”的图像转换为“myNameID123.nii”
    %d：插入系列描述（0008,103E）。例如，使用“myName％d”转换的回波平面图像会产生“myNameEPI”
    %s：插入系列（DICOM标签0020,0011）。例如，使用输出文件名“myName％s”将第二个序列转换为“myName2.nii”。
    %z：插入序列名称（0018,0024），因此使用“myName％z”转换的T1扫描可能会产生“myNameT1”。

    -o: 输出文件夹
    -z: 是否压缩

  • pydicom - Python’s pydicom package may be used to read and write DICOM files and convert them to NIfTI.

  • SPM– Statistical Parametric Mapping is a program that may be used to analyze images

  • ITK-Snap - is an image analysis toolkit that may be used to segment images

像 SPM 和 ITK-Snap 这样的一些工具更复杂，需要培训才能操作，而像 dcm2nii 和 dcm2niix 这样的工具是基于命令行的，操作简单。还要注意的是，这些应用程序中的大多数都希望图片在转换之前处于特定的格式和组织结构，比如特定的文件夹结构。

许多成像程序现在能够读取和写入NIfTI文件，因为这种格式在该领域已被广泛接受。由于得到许多计算机语言和库的支持，NIfTI数据可以在各种程序中轻松操作。例如，ITK-Snap、FSL、SPM和AFNI都能读取NIfTI文件（以.nii或.nii.gz结尾）。

1.2.1 NIfTI file format

  • the actual image data
  • header information

文件包含头文件（text data）和图片数据（binary data）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import nibabel as nib

# 读取NIfTI文件
nii_path = '/root/workspace/project/fmri_pre/data/009_MB1200_32C_shiyan1_niix/009_MB1200_32C_shiyan1Lulin-zhangyuxin-day_HB_202403131406249.nii.gz'
nii_img = nib.load(nii_path)

# 获取数据数组
data = nii_img.get_fdata()
print(data.shape)

# 查看nifti文件的仿射矩阵，得到坐标原点和体素大小
affine = nii_img.affine
# 提取体素大小
voxel_sizes = (affine[0, 0], affine[1, 1], affine[2, 2])

# 提取坐标原点（世界坐标系中）
origin = affine[:3, 3]

print("体素大小(mm):", voxel_sizes)
print("坐标原点(世界坐标):", origin)

# 获取NIfTI头文件
header = nii_img.header

# 打印头文件的内容
print(header)

1.3 fMRI技术原理

​ 在外界视觉刺激下，视觉皮层中神经细胞激活，血管耗氧量增加，脑血流速率增长，由于脑增长的速率远大于脑氧代谢率，导致神经激活处氧含量的供应远大于消耗，使得脱氧血红蛋白比例降低。由于血红蛋白是否含氧在磁性上稍有差别，含血红蛋白是抗磁性的，而脱血红蛋白是顺磁性的，在磁场中会引起自旋散相，导致磁共振信号在含氧的血液中显示高信号，而在脱氧的血液中显示低信号。神经活动引起的含氧血红蛋白比例升高，使得这种血氧水平依赖(Blood Oxygenation Level Dependent,BOLD) 信号增加。因此， BOLD 信号的改变一定程度上可以反映视觉皮层区域神经活动的情况，MRI是目前了解大脑功能最好的工具之一。相比神经活动，血流增加大约存在5至8秒的滞后，因此， MRI 的时间分辨率不高，导致fMRI实验数据的采集需要耗费较长的时间，制约了 fMRI 视觉信息解析研究的数据规模。

​ fMRI采集大脑的BOLD信号，通过三维重构等预处理过程，可以得到一系列三维大脑图像，每个三维图像都包含大量体素。体素是三维空间中成像的基本单元，与二维图像中的像素对应，每个体素的时间信号就是该体素对应视觉皮层空间区域内神经元群体活动引起的BOLD信号。当被试观看图像刺激时，相应的视觉皮层MRI体素会随时间发生规律性的变化。一个体素大约包含10个神经细胞，体素的BOLD信号是大量神经细胞的血液动力学响应的综合结果，反映了该区域内大量神经活动的综合情况。[^1]

1.4 图像刺激下的一般范式

​ 作为视觉刺激的图像数据集通常从公开自然图像数据库中选择，为了后续计算模型的训练和验证，通常把数据集划分为训练集和测试集。然后，通常按照随机选择的顺序，使用刺激呈现软件如 E-Prime 等，根据预先设定的呈现时间、呈现频率、休息间隙等播放给被试观看。通常一次图像呈现和数据采集被称作一次试验(trial)，如格兰特团队使用 200毫秒的频率和灰色背景交替呈现，持续 1 秒，然后呈现 3 秒的灰色背景，总共 4 秒时间，接着进行下一幅图像刺激的呈现试验。若干图像刺激的不断呈现通常作为一次扫描(run)，多次扫描组成一次实验，通常，整个数据库图像的呈现又分成若干次实验。另外，作为测试集的图像通常被重复多次呈现，用于提高 BOLD 信号的信噪比。

2. 数据预处理

​ 在将数据输入到模型之前，需要进行一系列预处理步骤，以提高数据质量并减少噪声的影响。对于fMRI数据，常见的预处理步骤包括：

  • Fuse and Conform : fMRI实验中通常需要多次扫描，包括功能扫描和结构扫描，有时还包括不同的条件和重复实验。不同的扫描可能有不同的分辨率和方向。为了能够将这些不同的数据集准确地叠加在一起，需要先将它们转换成统一的标准。这一步也是使数据准备进一步分析（如脑区分割、功能定位等）的基础。

  • Skull-stripping : 去除头壳

  • Spatial Normalization：如果尚未完成，将所有被试的脑图像变换到标准空间（如 MNI 空间）。

    标准空间，如蒙特利尔神经学研究所（Montreal Neurological Institute, MNI）空间，是神经影像学中用于数据分析和比较的一种共同参考框架。在神经科学研究中，研究者经常需要比较来自不同个体的大脑结构或功能图像。然而，由于每个人的大脑大小、形状和定位都有所不同，直接比较这些图像会非常困难。为了解决这个问题，科学家们开发了标准空间，即将个体的大脑图像变换到一个共同的参考框架内，以便进行有效的比较和分析。

    MNI空间是最广泛使用的标准空间之一，由蒙特利尔神经学研究所创建。它基于多个健康成年人大脑的平均形状和大小，提供了一个标准化的3D参考框架。MNI空间允许研究者将来自不同个体的大脑成像数据对齐到相同的三维坐标系统中，从而使比较和群体分析成为可能。

    将个体大脑数据变换到标准空间的过程通常涉及几个步骤，包括：

    1. 配准（Registration）：使用线性（刚体或仿射）和非线性变换将个体大脑图像对齐到标准空间模板。
    2. 归一化（Normalization）：调整图像的尺度，以匹配标准空间的大小和形状。
    3. 重采样（Resampling）：在变换到标准空间后，根据标准空间的分辨率调整图像的体素大小。

  • Brain Tissue Segmentation : 将脑成像数据中的不同组织类型（如灰质、白质和脑脊液）区分开的过程。这通常是通过分析MRI图像的强度值来实现的，因为不同的脑组织类型在MRI扫描中呈现不同的信号强度。

3. 获取特定脑区数据

1. 使用标准脑图谱获取掩码

​ 许多公开的脑图谱已经标注了不同的脑区，包括视觉区域。例如：

  • Automated Anatomical Labeling (AAL) 图谱：提供了大脑各个区域的标注，包括视觉皮层区域。
  • Harvard-Oxford Atlas：同样提供了详细的脑区划分。
  • Brodmann Areas：根据皮层细胞结构的不同，划分了大脑的多个区域，其中包括视觉处理区域。

使用这些图谱，你可以直接获取到预定义的视觉区域掩码。

2. 使用神经影像软件获取掩码

​ 软件如 FSL 、SPM 、或 FreeSurfer 等，都可以利用上述图谱来生成特定脑区的掩码。例如，使用 FSL 的步骤如下：

  • 打开FSLeyes。
  • 加载你的标准脑（例如，MNI152标准脑）。
  • 找到对应的图谱并应用到标准脑上。
  • 选择你感兴趣的视觉区域。
  • 导出该区域作为NIfTI格式的掩码文件。

3. 自定义绘制掩码

​ 如果公开的图谱不满足你的需求，你可以使用 MRIcron 、ITK-SNAP 或 FSLeyes 等工具手动绘制掩码。

  • 在这些软件中打开一个脑部MRI图像作为背景。
  • 手动在你感兴趣的视觉区域上绘制。
  • 保存你绘制的区域为掩码文件。

4. 提取视觉区域的数据

​ 一旦获得了视觉区域的掩码文件，就可以使用Python（例如，使用nibabel库）来读取fMRI数据和掩码文件，然后提取视觉区域的数据了。

1
2
3
4
5
6
7
8
9
10
11
12
13
import nibabel as nib
import numpy as np

# 加载fMRI数据和视觉区掩码
fmri_img = nib.load('fmri_data.nii.gz')
mask_img = nib.load('visual_mask.nii.gz')

# 获取数据
fmri_data = fmri_img.get_fdata()
mask_data = mask_img.get_fdata()

# 掩码应用到fMRI数据
visual_region_data = fmri_data[mask_data > 0]

4. 视觉解码模型

​ 基于 fMRI 的视觉信息解码模型能够根据视觉皮层 fMRI 体素响应预测相应自然图像刺激的内容，如图像场景语义类别、场景内容等，通常根据解析信息程度的不同，把视觉解码细分为视觉分类和视觉重构。

4.1 视觉分类模型

4.11简介

​ 视觉分类模型首先在 fMRI 数据训练集上训练视觉皮层 fMRI 体素响应到图像场景语义类别的映射，然后，能够根据测试集中新的视觉皮层 fMRI 体素响应预测被试所观看图像刺激的场景语义类别。 通常视觉皮层中体素数量较多，即视觉皮层区域 fMRI 体素响应向量维度较高，且 fMRI 数据集规模较小，直接进行模型训练容易导致过拟合。 因此，在视觉皮层大量体素响应输入模型之前，通常从视觉区域中的大量体素中 挑选一些重要的体素，降低用于后续视觉分类模型构建的输入维度 ，从而提高分类模型的准确性和鲁棒性。

样本少，特征高维，极易过拟合原因：由于特征数量众多，模型有能力学习训练数据中的复杂模式。然而，由于样本如此之少，其中许多模式很可能是噪声而不是与认知状态相关的真实信号。模型可能最终会记住这些噪声模式，这些模式无法很好地推广到新的、未见过的数据。

4.1.2 主要降维方法

  • 基于感兴趣区域

    基于感兴趣区域的体素降维方法根据视觉皮层的解剖学定位或者功能定位，选取某些感兴趣区域内的体素，输入视觉分类模型，这样可以大大降低输入的体素规模。通常选取大脑视觉皮层中与图像识别等相关的区域作为 ROI。另外，由于绝大多数的 fMRI 信号都是存在于灰质上，通常只保留灰质上的体素。

  • 基于主成分分析法

    不同与体素选择的方法，该方法利用主成分分析(Principal Component Analysis, PCA)对体素降维。通过线性变化把体素映射到另外一个低维空间，通过构造这个低维空间，使得体素投影到新的低维空间，且能够保持原有较多的信息，通过选择低维空间中投影数量值较大的前几个维度，保持原始体素对方差贡献最大的特征实现了对体素的降维。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset

# 假设数据已加载和预处理
# data的形状为[N, D, H, W, T]，其中N是样本数，D是深度，H是高度，W是宽度，T是时间点数
# labels是相应的标签，形状为[N]
# 为了简化，这里只考虑一个时间点的情况，所以我们使用data[:, :, :, :, 0]或者任何一个时间点

# 示例3D CNN模型
class Simple3DCNN(nn.Module):
    def __init__(self):
        super(Simple3DCNN, self).__init__()
        self.conv1 = nn.Conv3d(1, 32, kernel_size=3, stride=1, padding=1)  # 输入通道为1，输出通道为32
        self.pool = nn.MaxPool3d(kernel_size=2, stride=2, padding=0)
        self.conv2 = nn.Conv3d(32, 64, kernel_size=3, stride=1, padding=1)
        self.fc1 = nn.Linear(64 * 26 * 26 * 34, 128)  # 假设经过两次池化后的维度
        self.fc2 = nn.Linear(128, 2)  # 假设是一个二分类问题

    def forward(self, x):
        x = self.pool(F.relu(self.conv1(x)))
        x = self.pool(F.relu(self.conv2(x)))
        x = x.view(-1, 64 * 26 * 26 * 34)  # Flatten操作
        x = F.relu(self.fc1(x))
        x = self.fc2(x)
        return x

# 数据准备
# 假设你已有的data和labels
# data = ...  # 预处理后的fMRI数据，形状[N, 1, D, H, W]
# labels = ...  # 标签

# 转换为torch tensors
data_tensor = torch.Tensor(data)  # 需要将数据转换为torch张量
labels_tensor = torch.Tensor(labels).long()  # 标签也需要转换为torch张量，使用.long()转换为整数类型

# 创建数据加载器
dataset = TensorDataset(data_tensor, labels_tensor)
dataloader = DataLoader(dataset, batch_size=4, shuffle=True)

# 初始化模型、损失函数和优化器
model = Simple3DCNN()
criterion = nn.CrossEntropyLoss()
optimizer = optim.Adam(model.parameters(), lr=0.001)

# 训练模型
num_epochs = 10  # 训练周期数
for epoch in range(num_epochs):
    for i, (inputs, labels) in enumerate(dataloader):
        optimizer.zero_grad()
        outputs = model(inputs)
        loss = criterion(outputs, labels)
        loss.backward()
        optimizer.step()
        
    print(f'Epoch [{epoch+1}/{num_epochs}], Loss: {loss.item()}')

参考链接：

[2] DICOM 详解

[3] NifTi 介绍