1 % Batch SPM5 preprocessing (not all features obviouly)   R. Henson, Jan 06
   2 % assumes filenames code scan number
   3 % only explicitly sets non-default flag values
   4 
   5 owd = '/myhome';
   6 cd(owd)
   7 
   8 which spm
   9 spm('defaults', 'FMRI')
  10 global defaults
  11 
  12 subnum=[1:18];
  13 subnam = [050045 060001 060002 060003 060004 060005 060013 060014 060015 060016 060017 060018 060019 060023 060024 060025 060035 060036];
  14 
  15 sesdir = {'Ses1','Ses2','Deb'};
  16 
  17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  18 % Some flags to run control
  19 
  20 chkspikes = 1;
  21 realignflag = 1;
  22 %unwarpflag = 0;
  23 normviaT1 = 1;
  24 normviaT1seg = 1;
  25 
  26 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  27 
  28 for s = 1:length(subnum);
  29 
  30  sub = subnum(s) 
  31 
  32  if subnam(s)==050045 
  33 	nscan=[735];
  34  elseif subnam(s)==060004 
  35 	nscan=[735 735 150];
  36  else
  37 	nscan=[740 740 150];
  38  end
  39  nses = length(nscan);
  40  scans = cell(1,nses);
  41 
  42  swd = sprintf('%s/CBU%06d',owd,subnam(s))
  43  
  44  try, cd(swd), catch, error('No subject directory'), end
  45 
  46  for n = 1:nses
  47 
  48  	[files, flag] = spm_select('List',sesdir{n},'^fCBU.*');
  49 
  50   	if flag==0 | size(files,1) ~= nscan(n)
  51 		error(sprintf('Wrong no. of scans %d',size(files,1)))
  52 	else
  53 	  for f = 1:nscan(n)
  54 	   scans{n}=[scans{n}; fullfile(sesdir{n},files(f,:))];
  55 	  end
  56 	end
  57 
  58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  59 % Check for outlier slices (more than 5 stdevs)
  60 % note spm5spike is not part of SPM5
  61 % (if diff j is outlier, assumes scan j+1 at fault)
  62 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  63  
  64 	if chkspikes
  65 	 [aout,aout_scans,adout,adout_scans] = spm5spike(scans{n},6,1,1);
  66 	 eval(sprintf('print -djpeg -f1 spikes_sub%d_ses%d.jpg',sub,n)) 
  67 	 eval(sprintf('save spikes_sub%d_ses%d adout adout_scans',sub,n)) 
  68 	 outlier_nos{sub,n}=adout;
  69 	 outlier_scans{sub,n}=adout_scans;
  70 	end
  71  end	
  72 
  73 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  74 % Realignment (coregister & create mean only)
  75 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  76 
  77  if realignflag
  78   disp(sprintf('sub %d, realigning',sub))
  79 
  80   defaults.realign.estimate.quality = 1;
  81   defaults.realign.estimate.sep = 2;
  82   defaults.realign.estimate.wrap = [0 1 0];
  83   defaults.realign.estimate.interp = 7;
  84   defaults.realign.estimate.rtm = 1;
  85 
  86   spm_realign(scans,defaults.realign.estimate);
  87 
  88 % this is for reslicing mean only, do 'which',2, if want all
  89   defaults.realign.write.which = 0;
  90   defaults.realign.write.wrap = [0 1 0];
  91   defaults.realign.write.interp = 7;
  92   defaults.realign.write.mask = 1;
  93   defaults.realign.write.mean = 1;
  94 
  95   disp(sprintf('sub %d, creating mean...',sub))
  96 
  97   spm_reslice(strvcat(scans{:}),defaults.realign.write);
  98  end
  99 
 100 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 101 % Coreg EPIs and Struc
 102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 103 
 104  [struc,flag] = spm_select('List',fullfile(swd,'structurals'),'^sCBU.*');
 105  struc = fullfile(swd,'structurals',deblank(struc(1,:)));
 106 
 107  [meanf,flag] = spm_select('List',sesdir{1},'^mean.*');
 108  meanf = fullfile(sesdir{1},meanf);
 109 
 110 %%%%%%%%%%%%%%%%%%%%%%%%%
 111 % Coreg struc to mean func
 112 %
 113 % disp(sprintf('sub %d, coregistering structural to mean functional')) 
 114 %
 115 % defaults.coreg.estimate.cost_fun = 'nmi';
 116 %
 117 % x  = spm_coreg(spm_vol(meanf), spm_vol(struc), defaults.coreg.estimate);
 118 % M  = inv(spm_matrix(x));
 119 % MM = spm_get_space(deblank(struc));
 120 % spm_get_space(deblank(struc), M*MM);
 121 
 122 
 123 %%%%%%%%%%%%%%%%%%%%%%%%%
 124 % Coreg funcs to struc 
 125 %
 126 % disp(sprintf('sub %d, coregistering functionals to structural')) 
 127 %
 128 % defaults.coreg.estimate.cost_fun = 'nmi';
 129 %
 130 % x  = spm_coreg(spm_vol(struc), spm_vol(meanf), defaults.coreg.estimate);
 131 % M  = inv(spm_matrix(x));
 132 % MM = spm_get_space(deblank(meanf));
 133 % spm_get_space(deblank(meanf), M*MM);
 134 %
 135 % PO = strvcat(scans{:});
 136 % MM = zeros(4,4,size(PO,1));
 137 % for j=1:size(PO,1),
 138 %       MM(:,:,j) = spm_get_space(deblank(PO(j,:)));
 139 % end;
 140 % for j=1:size(PO,1),
 141 %       spm_get_space(deblank(PO(j,:)), M*MM(:,:,j));
 142 % end;
 143 
 144  if normviaT1
 145 
 146   if normviaT1seg
 147 
 148 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 149 % Determine normalisation params via Segmentation
 150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 151 
 152   disp(sprintf('sub %d, segmenting structural')) 
 153 
 154 %%%%%%%%%%%%%%%%%% 1st pass (to correct attenuation)
 155 
 156   estopts.regtype='';			% turn off affine
 157   out = spm_preproc(struc,estopts);
 158   
 159   [sn,isn]   = spm_prep2sn(out);
 160 
 161   % only write out attenuation corrected image
 162   writeopts.biascor = 1;
 163   writeopts.GM  = [0 0 0];
 164   writeopts.WM  = [0 0 0];
 165   writeopts.CSF = [0 0 0];
 166   writeopts.cleanup=0;
 167   spm_preproc_write(sn,writeopts);
 168 
 169   [pth,nam,ext,num] = spm_fileparts(struc);
 170   struc = fullfile(pth,['m' nam ext]);
 171 
 172 %%%%%%%%%%%%%%%%%% 2nd pass
 173 
 174   estopts.regtype='mni';		% turn on affine
 175   out = spm_preproc(struc,estopts);
 176 
 177   [sn,isn]   = spm_prep2sn(out);
 178 
 179   writeopts.biascor = 1;
 180   writeopts.GM  = [0 1 1];	% assume GM(2) means unmod
 181   writeopts.WM  = [0 0 0];
 182   writeopts.CSF = [0 0 0];
 183   spm_preproc_write(sn,writeopts);
 184 
 185   save(sprintf('%s_seg_sn.mat',spm_str_manip(struc,'sd')),'sn')
 186   save(sprintf('%s_seg_inv_sn.mat',spm_str_manip(struc,'sd')),'isn')
 187 
 188   [pth,nam,ext,num] = spm_fileparts(struc);
 189   struc = fullfile(pth,['m' nam ext]);
 190 
 191 %%%%%%%%%%%%%%%%%%%%%%%%%
 192 % Coreg funcs to atten corrected struc 
 193 % (REALLY WORTH RISK IF COREG FROM SCANNER OK (IF MINIMAL MOVEMENT?)
 194 
 195   disp(sprintf('sub %d, coregistering functionals to structural')) 
 196 
 197   defaults.coreg.estimate.cost_fun = 'nmi';
 198 
 199   x  = spm_coreg(spm_vol(struc), spm_vol(meanf), defaults.coreg.estimate);
 200   M  = inv(spm_matrix(x));
 201   MM = spm_get_space(deblank(meanf));
 202   spm_get_space(deblank(meanf), M*MM);
 203 
 204   PO = strvcat(scans{:});
 205   MM = zeros(4,4,size(PO,1));
 206   for j=1:size(PO,1),
 207        MM(:,:,j) = spm_get_space(deblank(PO(j,:)));
 208   end;
 209   for j=1:size(PO,1),
 210        spm_get_space(deblank(PO(j,:)), M*MM(:,:,j));
 211   end
 212 
 213   else
 214 
 215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 216 % Deteremine normalisation params direct from T1 template
 217 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 218 
 219    sn_name  = [spm_str_manip(struc,'sd') '_T1_sn.mat']
 220    sn = spm_normalise(fullfile(spm('Dir'),'templates','T1.nii'),spm_vol(struc),sn_name);
 221 
 222   end
 223 
 224  else
 225 
 226   sn_name  = [spm_str_manip(meanf,'sd') '_EPI_sn.mat']
 227   sn = spm_normalise(fullfile(spm('Dir'),'templates','EPI.nii'),spm_vol(meanf),sn_name);
 228 
 229  end
 230 
 231 
 232 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 233 % Normalise and Smooth functions
 234 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 235 
 236  defaults.normalise.write.interp = 7;
 237  defaults.normalise.write.wrap = [0 1 0];
 238  defaults.normalise.write.vox = [1 1 1];
 239 
 240  disp(sprintf('sub %d, writing normalised structural',sub))
 241 
 242  spm_write_sn(struc,sn,defaults.normalise.write);
 243 
 244  defaults.normalise.write.vox = [3 3 3];
 245 
 246  disp(sprintf('sub %d, writing normalised mean',sub))
 247 
 248  spm_write_sn(meanf,sn,defaults.normalise.write);
 249 
 250 % masking not necessary?
 251 % msk = spm_write_sn(strvcat(scans),params,defaults.normalise.write,'mask');
 252 % spm_write_sn(meanf,params,defaults.normalise.write,msk);
 253 
 254  disp(sprintf('sub %d, writing smoothed normalised EPIs',sub))
 255 
 256  for n=1:nses
 257    for m=1:length(scans{n})
 258      VO = spm_write_sn(scans{n}(m,:),sn,defaults.normalise.write);
 259 %     VO = spm_write_sn(scans{n}(m,:),sn,defaults.normalise.write,msk);
 260      [pth,nam,ext] = fileparts(scans{n}(m,:));
 261      spm_smooth(VO,fullfile(pth,['sw' nam ext]),[8 8 8]);
 262 %     [LASTMSG, LASTID] = lastwarn;
 263 %     warning('off',LASTID);
 264      warning off;
 265   end
 266  end
 267 
 268  warning on;
 269 end
 270 
 271 cd ..