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的开山鼻祖伦敦学院里流过来的，可以参考）

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

%-----------------------------------------------------------------------
% 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