基于Dpabi和spm12的脑脊液（csf）分割和提取笔记

一、前言

脑脊液（csf）一直被认为与新陈代谢有重要关联，其为许多神经科学研究提供重要价值，从fMRI图像中提取脑脊液信号可用于多种神经系统疾病的诊断。特别是自2019年Science上那篇著名的csf-BOLD文章发表后，大家都试图用csf阅读出更多的信息。

本人提取csf信号的思想都基于封面2019年这篇Science，不过因为扫描参数不同略有调整。我在学习csf-extract时还是花了一个月的时间，并且踩了不少坑，特此在这篇笔记里梳理与分享。

二、环境配置

系统：Windows 11
软件运行平台：Matlab2022a
所需工具包：

SPM12（下载安装教程）
DPABI（下载安装教程）

三、脑脊液提取过程

csf protocol

一、前言

二、环境配置

三、脑脊液提取过程

1.数据预处理

a.时间层校正（Slice Timing）

b.头动校正（Realign）

2.提取csf掩模

3.抓取csf信号

总结

1.数据预处理

DPARSF的基础教程在此不做赘述，如需要可以参考（基础教程）。

将功能像的数据都由dcm转化为.nii文件并去除了前面10个时间点后，就可以开始预处理了。

a.时间层校正（Slice Timing）
首先必须要进行的是Slice Timing，因为各扫描层激发时间不同Slice timing校正可以将所有切片的信号同步到一个参考时间点,通常是TR(重复时间)的中点。打开DPARSF后在红框所圈范围分别填上序列对应的Slice Number,Slice Order和参考切片（一般为中间）。然后点击右下角运行。
DPRASFslicetiming前后对比，会发现信号强度经过了细微调整。
b.头动校正（Realign）
【可选】头动校正（Realign）。有些用短回波成像的研究中并没有进行这一步，因为会改变csf附近的体素信号，但是我处理的是正常TR=2的序列。如果不进行头动校正影响较大，我还是选择了Realign，方法仅供参考。但无论是否应用头动校正，都需要DPARSF的头动参数.txt文件，让我们排除掉头动过大的被试。

【可选】滤波（Filter）。滤波参数我选择的是0.01-0.1HZ因为低频信号通常包含呼吸、心跳等生理噪声，高频信号可能包含仪器噪声等无关信息。
要不要做滤波？ -可以都试。很多作者也报告省略滤波后的结果。但我用没有bandpass的信号进行相关矩阵分析发现：如果不做带通滤波，被试间的信号都显著相关，说明噪音大；bandpass后被试的信号只和自己的序列相关，与他人无关，说明提升了信噪比。
如何做滤波？ -用DPARSF请手动检查！最好自己写脚本，报错少。DPARSF的滤波功能应该是使用y_bandpass(）函数完成的，会在某些情况下在一个时间点上生成全是0信号。如果用DPARSF自带的滤波，一定要在处理完后逐个时间点检查。

滤波前后

2.提取csf掩模

在前面的步骤中，我们将BOLD信号经过了时间校正（Slice Timing）、头动校正（Reaglign）、带通滤波（Filter）后已经得到了信噪比较高的BOLD图像。

对预处理后的4DBOLD进行时间维度求平均，得到单张3D均值图像，然后利用SPM12的分割工具对该均值图像进行组织分割，得到CSF的概率图。之后设置较高的概率阈值(0.9)创建CSF掩模，以确保提取的是高置信度的CSF区域。最后将这个掩模应用到原始的4D功能像上，即可得到目标区域内的CSF时间序列信号。这种方法的优势在于通过严格的阈值筛选，保证了提取信号的可靠性。

% 计算均值图像并用SPM分割得到CSF掩模
img = mean(func_4d.img, 4);
save_nii(img, 'mean.nii');
spm_segment('mean.nii'); 

% 提取CSF区域信号
mask = double(csf_prob.img)/256 >= 0.9;
csf_timeseries = func_4d.img .* repmat(mask,[1,1,1,size(func_4d.img,4)]);

3.抓取csf信号

脑脊液主要分布在脑室系统和蛛网膜下隙中。在大脑横断面水平切片上，从下到上看，最大的CSF池是位于脑干前方的桥前池，这是在底部切片的位置。这就是我们在提取CSF信号时主要关注最后非零切片的原因，在这个位置上的CSF信号最为纯净，受到周围灰质和白质组织的污染最少。

% 对每个时间点的3D体积提取CSF信号
for i = 1:size(img,4)
    % 获取当前时间点的3D体积
    img3d = img(:,:,:,i);
    
    % 找到含有信号的最后一个切片（对应脑底部CSF池的位置）
    nonZeroSlices = squeeze(any(any(img3d, 1), 2));
    lastNonZeroSlice = find(nonZeroSlices, 1);
    img2d = img3d(:,:,lastNonZeroSlice);
    csf(j,i) = mean(img2d(img2d>0));
end