% Batch SPM5 preprocessing (not all features obviouly)   R. Henson, Jan 06
% assumes filenames code scan number
% only explicitly sets non-default flag values

owd = '/myhome';
cd(owd)

which spm
spm('defaults', 'FMRI')
global defaults

subnum=[1:18];
subnam = [050045 060001 060002 060003 060004 060005 060013 060014 060015 060016 060017 060018 060019 060023 060024 060025 060035 060036];

sesdir = {'Ses1','Ses2','Deb'};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some flags to run control

chkspikes = 1;
realignflag = 1;
%unwarpflag = 0;
normviaT1 = 1;
normviaT1seg = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for s = 1:length(subnum);

 sub = subnum(s) 

 if subnam(s)==050045 
	nscan=[735];
 elseif subnam(s)==060004 
	nscan=[735 735 150];
 else
	nscan=[740 740 150];
 end
 nses = length(nscan);
 scans = cell(1,nses);

 swd = sprintf('%s/CBU%06d',owd,subnam(s))
 
 try, cd(swd), catch, error('No subject directory'), end

 for n = 1:nses

 	[files, flag] = spm_select('List',sesdir{n},'^fCBU.*');

  	if flag==0 | size(files,1) ~= nscan(n)
		error(sprintf('Wrong no. of scans %d',size(files,1)))
	else
	  for f = 1:nscan(n)
	   scans{n}=[scans{n}; fullfile(sesdir{n},files(f,:))];
	  end
	end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check for outlier slices (more than 5 stdevs)
% note spm5spike is not part of SPM5
% (if diff j is outlier, assumes scan j+1 at fault)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
	if chkspikes
	 [aout,aout_scans,adout,adout_scans] = spm5spike(scans{n},6,1,1);
	 eval(sprintf('print -djpeg -f1 spikes_sub%d_ses%d.jpg',sub,n)) 
	 eval(sprintf('save spikes_sub%d_ses%d adout adout_scans',sub,n)) 
	 outlier_nos{sub,n}=adout;
	 outlier_scans{sub,n}=adout_scans;
	end
 end	

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Realignment (coregister & create mean only)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 if realignflag
  disp(sprintf('sub %d, realigning',sub))

  defaults.realign.estimate.quality = 1;
  defaults.realign.estimate.sep = 2;
  defaults.realign.estimate.wrap = [0 1 0];
  defaults.realign.estimate.interp = 7;
  defaults.realign.estimate.rtm = 1;

  spm_realign(scans,defaults.realign.estimate);

% this is for reslicing mean only, do 'which',2, if want all
  defaults.realign.write.which = 0;
  defaults.realign.write.wrap = [0 1 0];
  defaults.realign.write.interp = 7;
  defaults.realign.write.mask = 1;
  defaults.realign.write.mean = 1;

  disp(sprintf('sub %d, creating mean...',sub))

  spm_reslice(strvcat(scans{:}),defaults.realign.write);
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Coreg EPIs and Struc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 [struc,flag] = spm_select('List',fullfile(swd,'structurals'),'^sCBU.*');
 struc = fullfile(swd,'structurals',deblank(struc(1,:)));

 [meanf,flag] = spm_select('List',sesdir{1},'^mean.*');
 meanf = fullfile(sesdir{1},meanf);

%%%%%%%%%%%%%%%%%%%%%%%%%
% Coreg struc to mean func
%
% disp(sprintf('sub %d, coregistering structural to mean functional')) 
%
% defaults.coreg.estimate.cost_fun = 'nmi';
%
% x  = spm_coreg(spm_vol(meanf), spm_vol(struc), defaults.coreg.estimate);
% M  = inv(spm_matrix(x));
% MM = spm_get_space(deblank(struc));
% spm_get_space(deblank(struc), M*MM);


%%%%%%%%%%%%%%%%%%%%%%%%%
% Coreg funcs to struc 
%
% disp(sprintf('sub %d, coregistering functionals to structural')) 
%
% defaults.coreg.estimate.cost_fun = 'nmi';
%
% x  = spm_coreg(spm_vol(struc), spm_vol(meanf), defaults.coreg.estimate);
% M  = inv(spm_matrix(x));
% MM = spm_get_space(deblank(meanf));
% spm_get_space(deblank(meanf), M*MM);
%
% PO = strvcat(scans{:});
% MM = zeros(4,4,size(PO,1));
% for j=1:size(PO,1),
%       MM(:,:,j) = spm_get_space(deblank(PO(j,:)));
% end;
% for j=1:size(PO,1),
%       spm_get_space(deblank(PO(j,:)), M*MM(:,:,j));
% end;

 if normviaT1

  if normviaT1seg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine normalisation params via Segmentation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  disp(sprintf('sub %d, segmenting structural')) 

%%%%%%%%%%%%%%%%%% 1st pass (to correct attenuation)

  estopts.regtype='';			% turn off affine
  out = spm_preproc(struc,estopts);
  
  [sn,isn]   = spm_prep2sn(out);

  % only write out attenuation corrected image
  writeopts.biascor = 1;
  writeopts.GM  = [0 0 0];
  writeopts.WM  = [0 0 0];
  writeopts.CSF = [0 0 0];
  writeopts.cleanup=0;
  spm_preproc_write(sn,writeopts);

  [pth,nam,ext,num] = spm_fileparts(struc);
  struc = fullfile(pth,['m' nam ext]);

%%%%%%%%%%%%%%%%%% 2nd pass

  estopts.regtype='mni';		% turn on affine
  out = spm_preproc(struc,estopts);

  [sn,isn]   = spm_prep2sn(out);

  writeopts.biascor = 1;
  writeopts.GM  = [0 1 1];	% assume GM(2) means unmod
  writeopts.WM  = [0 0 0];
  writeopts.CSF = [0 0 0];
  spm_preproc_write(sn,writeopts);

  save(sprintf('%s_seg_sn.mat',spm_str_manip(struc,'sd')),'sn')
  save(sprintf('%s_seg_inv_sn.mat',spm_str_manip(struc,'sd')),'isn')

  [pth,nam,ext,num] = spm_fileparts(struc);
  struc = fullfile(pth,['m' nam ext]);

%%%%%%%%%%%%%%%%%%%%%%%%%
% Coreg funcs to atten corrected struc 
% (REALLY WORTH RISK IF COREG FROM SCANNER OK (IF MINIMAL MOVEMENT?)

  disp(sprintf('sub %d, coregistering functionals to structural')) 

  defaults.coreg.estimate.cost_fun = 'nmi';

  x  = spm_coreg(spm_vol(struc), spm_vol(meanf), defaults.coreg.estimate);
  M  = inv(spm_matrix(x));
  MM = spm_get_space(deblank(meanf));
  spm_get_space(deblank(meanf), M*MM);

  PO = strvcat(scans{:});
  MM = zeros(4,4,size(PO,1));
  for j=1:size(PO,1),
       MM(:,:,j) = spm_get_space(deblank(PO(j,:)));
  end;
  for j=1:size(PO,1),
       spm_get_space(deblank(PO(j,:)), M*MM(:,:,j));
  end

  else

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Deteremine normalisation params direct from T1 template
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   sn_name  = [spm_str_manip(struc,'sd') '_T1_sn.mat']
   sn = spm_normalise(fullfile(spm('Dir'),'templates','T1.nii'),spm_vol(struc),sn_name);

  end

 else

  sn_name  = [spm_str_manip(meanf,'sd') '_EPI_sn.mat']
  sn = spm_normalise(fullfile(spm('Dir'),'templates','EPI.nii'),spm_vol(meanf),sn_name);

 end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Normalise and Smooth functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 defaults.normalise.write.interp = 7;
 defaults.normalise.write.wrap = [0 1 0];
 defaults.normalise.write.vox = [1 1 1];

 disp(sprintf('sub %d, writing normalised structural',sub))

 spm_write_sn(struc,sn,defaults.normalise.write);

 defaults.normalise.write.vox = [3 3 3];

 disp(sprintf('sub %d, writing normalised mean',sub))

 spm_write_sn(meanf,sn,defaults.normalise.write);

% masking not necessary?
% msk = spm_write_sn(strvcat(scans),params,defaults.normalise.write,'mask');
% spm_write_sn(meanf,params,defaults.normalise.write,msk);

 disp(sprintf('sub %d, writing smoothed normalised EPIs',sub))

 for n=1:nses
   for m=1:length(scans{n})
     VO = spm_write_sn(scans{n}(m,:),sn,defaults.normalise.write);
%     VO = spm_write_sn(scans{n}(m,:),sn,defaults.normalise.write,msk);
     [pth,nam,ext] = fileparts(scans{n}(m,:));
     spm_smooth(VO,fullfile(pth,['sw' nam ext]),[8 8 8]);
%     [LASTMSG, LASTID] = lastwarn;
%     warning('off',LASTID);
     warning off;
  end
 end

 warning on;
end

cd ..