SPM12之fMRI批量预处理——NII文件处理
主要流程：Realignment、Slice Timing Correction、Coregistration、Segmentation、Normalization
代码地址：https://github.com/zhoushusheng/fMRI_preprocessing
参数和步骤设置参考：
1、上海外国网课https://www.bilibili.com/video/BV1ye4y1m72P/spm_id_from=333.880.my_history.page.click&vd_source=9b75cfd2936571e6bc329144e79be03e（这里学处理步骤的设置，视频里说预处理有两种方式：1、在Coregister里，依照Realignment生成的mean文件，选择自身的MRI进行coregistration。2、选择的不是自己的MRI，好像是用统一模板，还是说直接跳过了，normalization里就不用Coregistration的文件，读者这里自己再考证）
2、Andy's brain book：https://andysbrainbook.readthedocs.io/en/latest/SPM/SPM_Short_Course/SPM_04_Preprocessing.html（这个很权威，一些的参数设置是这里的）
3. slice order : 最难确定的就是SLICE ORDER了。查阅了很多资料，一般包含在你下载的DICOM文件里，但我没发现。链接里说（https://rfmri.org/content/slice-order-adni-rsfmri-data）：或者需要机器型号才能确定。实在不知道，就得问核磁共振的操作员了（我上哪找，这不为难人吗）。上海外国网课https://www.bilibili.com/video/BV1ye4y1m72P/spm_id_from=333.880.my_history.page.click&vd_source=9b75cfd2936571e6bc329144e79be03e 里说slice order 要么是从上往下，要么从下往上，第一层一般是在底部，扫描有可能奇数偶数交叉。参数可能都试过了，各种参数对最后出来的结果影像不大（通过生成皮尔逊ROI相关矩阵判断），我初步推断因为SLICE ORDER是Slice Timing Correction的操作，且核磁共振扫描是一层层的扫描，所以呈相会有时间差，但都是几十或者上千毫秒（此话出处：Andy's brain book：https://andysbrainbook.readthedocs.io/en/latest/SPM/SPM_Short_Course/SPM_04_Preprocessing.html），最后可以推断slice order 是错的也无伤大雅.。这是此问题的回答锦集，虽然没有直接给出答案：https://blog.csdn.net/Junbao_Zhang/article/details/85245922。
个人的批量代码的参数生成教程：https://blog.csdn.net/Iris_bysshqx17/article/details/100188483?spm=1001.2101.3001.6661.1&utm_medium=distribute.pc_relevant_t0.none-task-blog-2~default~BlogCommendFromBaidu~Rate-1-100188483-blog-125808030.pc_relevant_vip_default&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-2~default~BlogCommendFromBaidu~Rate-1-100188483-blog-125808030.pc_relevant_vip_default&utm_relevant_index=1（没有用这篇博客的参数设置和步骤，用的是他生成代码的过程，代码是那抄的，有错误，本文已改）
参考资料：
北师大的PPT:https://wenku.baidu.com/view/ea412f804afe04a1b071de7f.html?_wkts_=1721898330514（北师大的fMRI研究是国内领先的，他的PPT 也是SPM的开山鼻祖伦敦学院里流过来的，可以参考）
%-----------------------------------------------------------------------
% Job saved on 02-Sep-2019 18:21:16 by cfg_util (rev $Rev: 6942 $)
% spm SPM - SPM12 (7219)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
%%
% 这个m文件用来运行spm的matlabbatch
% 该文件中的文件位置是基于linux系统的，如果是windows系统需要修改/为\
clear
% spm_path = '~\spm12'; %设置spm12的位置（这里填的是安装=spm12的绝对位置） 我这里已经安装好了
% addpath(spm_path);
spm('defaults', 'fmri');
spm_jobman('initcfg');
%%
subjectsdir = {'B:\fMRI\preprocess_for_batch\CN\'};% 这里是data文件夹的绝对位置
subjects = {'sbj01','sbj02','sbj03','sbj04','sbj05','sbj06','sbj07','sbj08','sbj09','sbj10','sbj11','sbj12','sbj13','sbj14','sbj15','sbj16','sbj17','sbj18','sbj19','sbj20','sbj21','sbj22','sbj23','sbj24','sbj25'...
,'sbj26','sbj27','sbj28','sbj29','sbj30','sbj31','sbj32','sbj33','sbj34','sbj35','sbj36'}; % 单个或多个被试的文件夹
funcdir = fullfile('fMRI'); % 第一个Session的文件夹（由于preprocessing.mat中只添加了一个Session，所以这里也只使用一个Session。）
% F = fullfile(FOLDERNAME1, FOLDERNAME2, ..., FILENAME) builds a full
% file specification F from the folders and file name specified. Input
% arguments FOLDERNAME1, FOLDERNAME2, etc. and FILENAME can be strings,
% character vectors, or cell arrays of character vectors. Non-scalar
% strings and cell arrays of character vectors must all be the same size.
% FULLFILE collapses inner repeated file separators unless they appear at
% the beginning of the full file specification. FULLFILE also collapses
% relative folders indicated by the dot symbol, unless they appear at
% the end of the full file specification. Relative folders indicated
% by the double-dot symbol are not collapsed.
%
% funcdir2 = fullfile('Session2'); % 第二个Session的文件夹
anatdir = fullfile('MRI'); % 结构像的文件夹
%%
% 设置基本的参数，在这里设置方便后期修改
% head = spm_vol(sessionfiles{1});
% dimension = head.dim;
% slice_number = dimension(1,3);
% slice_number = 35;
% TR = 3;
% TA = TR-TR\slice_number;
% slice_order = [1:2:slice_number 2:2:slice_number];
% refslice1 = 35;
% voxel_size= [3.75 3.75 4];
% smoothing_kernel = [6 6 6];
%%
%每一个被试都建立一个工作流程
nsubj = length(subjects); % 被试的数量
jobs = cell(nsubj,1); % job的数量，每个被试一个job
nrun=1; % session的数量
%ntimepoint=190; %删除了前面的10个时间点后剩余的时间点
%%
for csubj = 24:36
subjdir = {spm_select('CPath',subjects{csubj}, subjectsdir)};%Cpath 的c = canonical
% 结构像文件夹
adir = spm_select('CPath', anatdir, subjdir); %选择当前被试（csubj）的结构像文件夹
anatfile = strcat(spm_select('FPList', adir, '\\*.nii'),',1'); %在这个文件夹中.nii结尾的文件，即结构像文件 % As above, but return files with full paths (i.e. prefixes 'direc' to each)
%如果没有结构像文件就报错
if isequal(anatfile, '')
warning(sprintf('No anat file found for %s', ...
subjects{csubj}))
return
end
% Session1文件夹
fdir = {spm_select('CPath', funcdir, subjdir)}; %选择当前被试（csubj）的Session1文件夹
ffiles = spm_select('List', fdir, '\\*.nii'); %选择这个文件夹中所有.nii结尾的文件，即Session1的所有原始功能像文件
%如果没有功能像文件就报错
nimage = size(ffiles,1);
if nimage == 0
warning(sprintf('No functional file found for %s', subjects{csubj}))
return
end
cffiles = cellstr(ffiles);
ntimepoint = length(cffiles);
funcfiles = cell(1, nrun);
sessionfiles = cell(nrun,ntimepoint);
for i = 1:nrun
for j = 1:ntimepoint
sessionfiles{i,j}=strcat(fdir{1},'\',cffiles{j}); % 在这里在每一个功能像文件的绝对路径后面添加,1
end
funcfiles{1,i} = {sessionfiles{i,:}}'; % 如果有多个Sessions，funcfiles中就会有多个sessionfiles
end
clear matlabbatch
length(cffiles);
head = spm_vol(sessionfiles{1});
dimension = head.dim;
descrip = head.descrip;
slice_number = dimension(1,3);
descrip={descrip};%char To cell
TR = regexp(descrip,'(?<=TR=)\d+|(?<=TR=)\d+\.\d+','match','once');
TR = str2double(TR);
TR = TR / 1000;
%如果没有查到TR就报错
if isequal(TR, '')
warning(sprintf('No TR parameter found for %s', ...
subjects{csubj}))
return
end
TA = TR-TR/slice_number;
%这个就得自己看了
% slice_order = [1:2:slice_number 2:2:slice_number];
% slice_order = [slice_number:-1:1];
slice_order = [slice_number:-1:1];
% preprocessing job
display 'Creating preprocessing job' %在matlab的命令行窗口显示Creating preprocessing job
matlabbatch{1}.spm.spatial.realign.estwrite.data = {funcfiles{1,1}(:,1)};%这里不加1 也行 ，就是上面查找的时候
%%
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.quality = 0.9;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.sep = 4;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.fwhm = 5;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.rtm = 1;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.interp = 2;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.wrap = [0 0 0];
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.weight = '';
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.which = [2 1];
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.interp = 4;
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.wrap = [0 0 0];
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.mask = 1;
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.prefix = 'r';
matlabbatch{2}.spm.temporal.st.scans{1}(1) = cfg_dep('Realign: Estimate & Reslice: Realigned Images (Sess 1)', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','sess', '()',{1}, '.','cfiles'));
matlabbatch{2}.spm.temporal.st.nslices = slice_number;
matlabbatch{2}.spm.temporal.st.tr = TR;
matlabbatch{2}.spm.temporal.st.ta = TA;
matlabbatch{2}.spm.temporal.st.so = slice_order;
matlabbatch{2}.spm.temporal.st.refslice = slice_number/2;
matlabbatch{2}.spm.temporal.st.prefix = 'a';
matlabbatch{3}.spm.spatial.coreg.estimate.ref(1) = cfg_dep('Realign: Estimate & Reslice: Mean Image', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','rmean'));
matlabbatch{3}.spm.spatial.coreg.estimate.source = {anatfile};%ananfile 字符串末尾得加一个1，细看对比stupid_batch_nii_third_formal_compilable
matlabbatch{3}.spm.spatial.coreg.estimate.other = {''};
matlabbatch{3}.spm.spatial.coreg.estimate.eoptions.cost_fun = 'nmi';
matlabbatch{3}.spm.spatial.coreg.estimate.eoptions.sep = [4 2];
matlabbatch{3}.spm.spatial.coreg.estimate.eoptions.tol = [0.02 0.02 0.02 0.001 0.001 0.001 0.01 0.01 0.01 0.001 0.001 0.001];
matlabbatch{3}.spm.spatial.coreg.estimate.eoptions.fwhm = [7 7];
matlabbatch{4}.spm.spatial.preproc.channel.vols(1) = cfg_dep('Coregister: Estimate: Coregistered Images', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','cfiles'));
matlabbatch{4}.spm.spatial.preproc.channel.biasreg = 0.001;
matlabbatch{4}.spm.spatial.preproc.channel.biasfwhm = 60;
matlabbatch{4}.spm.spatial.preproc.channel.write = [0 1];
matlabbatch{4}.spm.spatial.preproc.tissue(1).tpm = {'C:\Program Files\MATLAB\R2021b\toolbox\spm12\tpm\TPM.nii,1'};
matlabbatch{4}.spm.spatial.preproc.tissue(1).ngaus = 1;
matlabbatch{4}.spm.spatial.preproc.tissue(1).native = [1 0];
matlabbatch{4}.spm.spatial.preproc.tissue(1).warped = [0 0];
matlabbatch{4}.spm.spatial.preproc.tissue(2).tpm = {'C:\Program Files\MATLAB\R2021b\toolbox\spm12\tpm\TPM.nii,2'};
matlabbatch{4}.spm.spatial.preproc.tissue(2).ngaus = 1;
matlabbatch{4}.spm.spatial.preproc.tissue(2).native = [1 0];
matlabbatch{4}.spm.spatial.preproc.tissue(2).warped = [0 0];
matlabbatch{4}.spm.spatial.preproc.tissue(3).tpm = {'C:\Program Files\MATLAB\R2021b\toolbox\spm12\tpm\TPM.nii,3'};
matlabbatch{4}.spm.spatial.preproc.tissue(3).ngaus = 2;
matlabbatch{4}.spm.spatial.preproc.tissue(3).native = [1 0];
matlabbatch{4}.spm.spatial.preproc.tissue(3).warped = [0 0];
matlabbatch{4}.spm.spatial.preproc.tissue(4).tpm = {'C:\Program Files\MATLAB\R2021b\toolbox\spm12\tpm\TPM.nii,4'};
matlabbatch{4}.spm.spatial.preproc.tissue(4).ngaus = 3;
matlabbatch{4}.spm.spatial.preproc.tissue(4).native = [1 0];
matlabbatch{4}.spm.spatial.preproc.tissue(4).warped = [0 0];
matlabbatch{4}.spm.spatial.preproc.tissue(5).tpm = {'C:\Program Files\MATLAB\R2021b\toolbox\spm12\tpm\TPM.nii,5'};
matlabbatch{4}.spm.spatial.preproc.tissue(5).ngaus = 4;
matlabbatch{4}.spm.spatial.preproc.tissue(5).native = [1 0];
matlabbatch{4}.spm.spatial.preproc.tissue(5).warped = [0 0];
matlabbatch{4}.spm.spatial.preproc.tissue(6).tpm = {'C:\Program Files\MATLAB\R2021b\toolbox\spm12\tpm\TPM.nii,6'};
matlabbatch{4}.spm.spatial.preproc.tissue(6).ngaus = 2;
matlabbatch{4}.spm.spatial.preproc.tissue(6).native = [0 0];
matlabbatch{4}.spm.spatial.preproc.tissue(6).warped = [0 0];
matlabbatch{4}.spm.spatial.preproc.warp.mrf = 1;
matlabbatch{4}.spm.spatial.preproc.warp.cleanup = 1;
matlabbatch{4}.spm.spatial.preproc.warp.reg = [0 0.001 0.5 0.05 0.2];
matlabbatch{4}.spm.spatial.preproc.warp.affreg = 'mni';
matlabbatch{4}.spm.spatial.preproc.warp.fwhm = 0;
matlabbatch{4}.spm.spatial.preproc.warp.samp = 3;
matlabbatch{4}.spm.spatial.preproc.warp.write = [0 1];
matlabbatch{4}.spm.spatial.preproc.warp.vox = NaN;
matlabbatch{4}.spm.spatial.preproc.warp.bb = [NaN NaN NaN
NaN NaN NaN];
matlabbatch{5}.spm.spatial.normalise.write.subj.def(1) = cfg_dep('Segment: Forward Deformations', substruct('.','val', '{}',{4}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','fordef', '()',{':'}));
matlabbatch{5}.spm.spatial.normalise.write.subj.resample(1) = cfg_dep('Realign: Estimate & Reslice: Resliced Images (Sess 1)', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','sess', '()',{1}, '.','rfiles'));
matlabbatch{5}.spm.spatial.normalise.write.woptions.bb = [-78 -112 -70
78 76 85];
matlabbatch{5}.spm.spatial.normalise.write.woptions.vox = [2 2 2];
matlabbatch{5}.spm.spatial.normalise.write.woptions.interp = 4;
matlabbatch{5}.spm.spatial.normalise.write.woptions.prefix = 'w';
matlabbatch{6}.spm.spatial.smooth.data(1) = cfg_dep('Normalise: Write: Normalised Images (Subj 1)', substruct('.','val', '{}',{5}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('()',{1}, '.','files'));
matlabbatch{6}.spm.spatial.smooth.fwhm = [8 8 8];
matlabbatch{6}.spm.spatial.smooth.dtype = 0;
matlabbatch{6}.spm.spatial.smooth.im = 0;
matlabbatch{6}.spm.spatial.smooth.prefix = 's';
% 保存matlabbatch
matfile = sprintf('preprocess_nii_%s.mat', subjects{csubj});
% save(matfile,'matlabbatch');
jobs{csubj} = matfile; %jobs这个变量中存储了所有被试的matlabbatch
['Start preprocessing job' num2str(csubj)]
spm_jobman('run', matlabbatch);
end
%%
% 执行预处理的工作
% for csubj = 15:20
% display 'Start preprocessing job' %在matlab的命令行窗口显示Start preprocessing job
% load(jobs{csubj});
% spm_jobman('run', matlabbatch);
% end
本文来自博客园，作者：z_s_s，转载请注明原文链接：https://www.cnblogs.com/zhoushusheng/p/18322099
浙公网安备 33010602011771号