% Batch script for preprocessing and statistics for single subject 
% event-related fMRI example dataset			R. Henson, 20/5/03

spm('defaults','FMRI')
global defaults
global Ufp
Ufp = 0.001;

nsub = 1;
nses = 1;

cwd = '/windows/D/SPMCourseData/SingleSub/SPM2'		% CHANGE TO YOUR PATH
cd(cwd)

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


for sub=1:nsub

 P = cell(1,nses);
 for ses=1:nses
  dir{ses} = sprintf('Data');		% directories can be generalised for >1 ses/sub
  P{ses}   = spm_get('Files',dir{ses},'sM*.img');
 end

 cd(dir{1});				% where spm2.ps is output

 V = spm_vol(P);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Realignment (coregister only)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 disp(sprintf('sub %d, coregistering...',sub))

 spm_realign(V);

% If want HIGHER quality (but slower!))
% FlagsC  = struct('quality',1,'fwhm',5,'rtm',1,'wrap',[0 1 0],'interp',7);
% spm_realign(V,FlagsC);

% eval('!ghostview spm2.ps')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unwarp (and reslice)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 disp(sprintf('sub %d, unwarping...',sub))

 uwe_flags = struct(  'order',       defaults.unwarp.estimate.basfcn,...
	              'sfield',      [],...
	              'regorder',    defaults.unwarp.estimate.regorder,...
	              'lambda',      defaults.unwarp.estimate.regwgt,...
	              'jm',          defaults.unwarp.estimate.jm,...
	              'fot',         [4 5],...		% pitch and roll
	              'sot',         [],...		% no second-order 
	              'fwhm',        defaults.unwarp.estimate.fwhm,...
	              'rem',         defaults.unwarp.estimate.rem,...
	              'noi',         defaults.unwarp.estimate.noi,...
	              'exp_round',   defaults.unwarp.estimate.expround);
	
 uwr_flags = struct(  'wrap',        [0 1 0],...
	              'interp',      7,...
	              'mask',        1,...
	              'which',       2,...
	              'mean',        1);

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Estimate unwarping parameters
   
 for j=1:length(P)
	tmpP 	        = spm_vol(P{j}(1,:));
	uwe_flags.M 	= tmpP.mat;
	
	ds 	= spm_uw_estimate(P{j},uwe_flags);
	ads(j) 	= ds;
	[path,name,ext,ver] = fileparts(P{j}(1,:));
	pefile = fullfile(path,[name '_uw.mat']);
	save(pefile,'ds');
 end
	
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Write unwarped images (more efficient if spm_uw_apply returned VO)

 spm_uw_apply(ads,uwr_flags);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Reload unwarped images

 cd(cwd);
 P = cell(1,nses);
 for ses=1:nses
  dir{ses} = sprintf('Data');		% directories can be generalised for >1 ses/sub
  P{ses}   = spm_get('Files',dir{ses},'usM*.img');
 end

 cd(dir{1});				% where spm2.ps is output

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Memory map the mean image (not affected by slice-timing below)

 meanf   = [spm_str_manip(P{1}(1,:),'h') '/mean',...
			      spm_str_manip(P{1}(1,:),'t')];
 Vm      = spm_vol(meanf);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Slice-timing correction (optional; could be before spatial realignment if interleaved)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 TR = 2;
 num_slices = 24;
 ref_slice  = 12;			    % middle slice in time
 slice_time = TR /num_slices;
 slice_gap  = slice_time;		    % continuous scanning
 slice_ord  = [num_slices:-1:1];			% descending acquisition [num_slices:-1:1]

 spm_slice_timing(P, slice_ord, ref_slice, [slice_time slice_gap]);

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Reload slice-time corrected images

 cd(cwd);
 P = cell(1,nses);
 for ses=1:nses
  dir{ses} = sprintf('Data');		% directories can be generalised for >1 ses/sub
  P{ses}   = spm_get('Files',dir{ses},'ausM*.img');
 end

 V       = spm_vol(P);
 V       = cat(1,V{:});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Normalization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Estimate unwarping parameters
 
 matname  = [spm_str_manip(Vm.fname,'sd') '_sn.mat'];
 VG       = fullfile(spm('Dir'),'templates','EPI.mnc');
 params   = spm_normalise(VG,Vm,matname,'','',defaults.normalise.estimate);

 msk      = spm_write_sn(V,params,defaults.normalise.write,'mask');

 disp(sprintf('sub %d, determining normalisation parameters...',sub))
 spm_write_sn(Vm,params,defaults.normalise.write,msk);		% write nmean

% eval('!ghostview spm2.ps')

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Write normalised

 disp(sprintf('sub %d, writing normalised',sub))
 for i=1:length(V),
        VO(i) = spm_write_sn(V(i),params,defaults.normalise.write,msk);
 end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Smoothing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 disp(sprintf('sub %d, smoothing',sub))
 for i=1:length(VO),
      [pth,nam,ext] = fileparts(V(i).fname);
      fname         = fullfile(pth,['sw' nam ext]);
      spm_smooth(VO(i),fname,8);
 end;

 clear VO

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Categorical Stats model (see README) for one session/subject
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cd(cwd)
load sots

swd = [cwd '/CatStats'];
eval(sprintf('!mkdir %s',swd));
cd(swd);


% number of scans and sessions
%---------------------------------------------------------------------------
SPM.nscan          = ones(1,nses)*351;

% basis functions and timing parameters
%---------------------------------------------------------------------------

SPM.xBF.name       = 'hrf (with time and dispersion derivatives)';
SPM.xBF.order      = 3;                
SPM.xBF.length     = 32;                % length in seconds
SPM.xBF.T          = num_slices;       	% number of time bins per scan
SPM.xBF.T0         = ref_slice;		% middle slice/timebin
SPM.xBF.UNITS      = 'scans';           % OPTIONS: 'scans'|'secs' for onsets
SPM.xBF.Volterra   = 1;                 % OPTIONS: 1|2 = order of convolution

% Trial specification: Onsets, duration (UNITS) and parameters for modulation
%---------------------------------------------------------------------------

ncon = 4;
cnam = {'N1','N2','F1','F2'};

for ses=1:nses
	for c=1:ncon
 		SPM.Sess(ses).U(c).name      = {cnam{c}};
 		SPM.Sess(ses).U(c).ons       = sots{c};
 		SPM.Sess(ses).U(c).dur       = 0;
 		SPM.Sess(ses).U(c).P(1).name = 'none';
	end
end

% design (user specified covariates, eg movement parameters)
%---------------------------------------------------------------------------
% (note inclusion of movement parameters makes unwarping somewhat redundant,
% but done here to illustrate facilities of SPM)

rnam = {'X','Y','Z','x','y','z'};
for ses=1:nses
	dir{ses} = '../Data';		% directories can be generalised for >1 ses/sub
	fn = spm_get('Files',dir{ses},'rp_*.txt');
	[r1,r2,r3,r4,r5,r6] = textread(fn,'%f%f%f%f%f%f');
	SPM.Sess(ses).C.C = [r1 r2 r3 r4 r5 r6];
	SPM.Sess(ses).C.name = rnam;
end

% global normalization: OPTIONS:'Scaling'|'None'
%---------------------------------------------------------------------------
SPM.xGX.iGXcalc    = 'None';

% low frequency confound: high-pass cutoff (secs) [Inf = no filtering]
%---------------------------------------------------------------------------
for ses=1:nses
	SPM.xX.K(ses).HParam = 128;
end

% intrinsic autocorrelations: OPTIONS: 'none'|'AR(1) + w'
%-----------------------------------------------------------------------
SPM.xVi.form       = 'AR(1)';


% specify data: matrix of filenames and TR
%===========================================================================

for ses=1:nses
	dir{ses} = '../Data';		% directories can be generalised for >1 ses/sub
	tmp{ses} = spm_get('files',dir{ses},'swau*.img');
end
SPM.xY.P = cat(1,tmp{:});

SPM.xY.RT          = TR;


% Configure design matrix
%===========================================================================
SPM = spm_fmri_spm_ui(SPM);


% Estimate parameters
%===========================================================================
SPM = spm_spm(SPM);


% Add extra contrasts
%===========================================================================

ocon = length(SPM.xCon);			% default contrasts

eb = eye(SPM.xBF.order);			% basis functions


cnam{1} = 'Canonical HRF: Faces > Baseline';
cwgt{1} = [kron([1 1 1 1],eb(1,:)) zeros(1,6)];		% zeros = movement params
ctyp{1} = 'T';

cnam{2} = 'Effects of REAL interest';
cwgt{2} = [eye(4*SPM.xBF.order) zeros(4*SPM.xBF.order,6)];			
ctyp{2} = 'F';

cnam{3} = 'Canonical HRF: F vs N';
cwgt{3} = [kron([-1 -1 1 1],eb(1,:)) zeros(1,6)];		
ctyp{3} = 'T';

cnam{4} = 'Canonical + Derivatives: 1 vs 2';
cwgt{4} = [kron([1 -1 1 -1],eb) zeros(SPM.xBF.order,6)];			
ctyp{4} = 'F';

cnam{5} = 'Movement parameters';
cwgt{5} = [zeros(6,4*SPM.xBF.order) eye(6)];
ctyp{5} = 'F';

for c = 1:length(cnam)

	cw = [cwgt{c} zeros(size(cwgt{c},1),1)]';	% pad with zero for constant

	SPM.xCon(c+ocon)     = spm_FcUtil('Set',cnam{c},ctyp{c},'c',cw,SPM.xX.xKXs);
end


% and evaluate
%---------------------------------------------------------------------------
spm_contrasts(SPM);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now press RESULTS (and see README)!


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Parametric Stats model (see README) for one session/subject
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cd(cwd)
load lags

swd = [cwd '/ParStats'];
eval(sprintf('!mkdir %s',swd));
cd(swd);


% number of scans and sessions
%---------------------------------------------------------------------------
SPM.nscan          = ones(1,nses)*351;

% basis functions and timing parameters
%---------------------------------------------------------------------------

SPM.xBF.name       = 'hrf (with time and dispersion derivatives)';
SPM.xBF.order      = 3;                
SPM.xBF.length     = 32;                % length in seconds
SPM.xBF.T          = num_slices;       	% number of time bins per scan
SPM.xBF.T0         = ref_slice;		% middle slice/timebin
SPM.xBF.UNITS      = 'scans';           % OPTIONS: 'scans'|'secs' for onsets
SPM.xBF.Volterra   = 1;                 % OPTIONS: 1|2 = order of convolution

% Trial specification: Onsets, duration (UNITS) and parameters for modulation
%---------------------------------------------------------------------------

for ses=1:nses
	SPM.Sess(ses).U(1).name      = {'Face'};
	SPM.Sess(ses).U(1).ons       = sots';
	SPM.Sess(ses).U(1).dur       = 0;

	SPM.Sess(ses).U(1).P(1).name = 'Lag';
	SPM.Sess(ses).U(1).P(1).P    = lag';
	SPM.Sess(ses).U(1).P(1).h    = 1;
	SPM.Sess(ses).U(1).P(1).i    = [1 2];

	SPM.Sess(ses).U(1).P(2).name = 'Fame';
	SPM.Sess(ses).U(1).P(2).P    = fam';
	SPM.Sess(ses).U(1).P(2).h    = 1;
	SPM.Sess(ses).U(1).P(2).i    = [1 3];

	SPM.Sess(ses).U(1).P(3).name = 'FameXLag';
	SPM.Sess(ses).U(1).P(3).P    = famlag';
	SPM.Sess(ses).U(1).P(3).h    = 1;
	SPM.Sess(ses).U(1).P(3).i    = [1 4];
end

% design (user specified covariates, eg movement parameters)
%---------------------------------------------------------------------------
% (note inclusion of movement parameters makes unwarping somewhat redundant,
% but done here to illustrate facilities of SPM)

rnam = {'X','Y','Z','x','y','z'};
for ses=1:nses
	dir{ses} = '../Data';		% directories can be generalised for >1 ses/sub
	fn = spm_get('Files',dir{ses},'rp_*.txt');
	[r1,r2,r3,r4,r5,r6] = textread(fn,'%f%f%f%f%f%f');
	SPM.Sess(ses).C.C = [r1 r2 r3 r4 r5 r6];
	SPM.Sess(ses).C.name = rnam;
end

% global normalization: OPTIONS:'Scaling'|'None'
%---------------------------------------------------------------------------
SPM.xGX.iGXcalc    = 'None';

% low frequency confound: high-pass cutoff (secs) [Inf = no filtering]
%---------------------------------------------------------------------------
for ses=1:nses
	SPM.xX.K(ses).HParam = 128;
end

% intrinsic autocorrelations: OPTIONS: 'none'|'AR(1) + w'
%-----------------------------------------------------------------------
SPM.xVi.form       = 'AR(1)';


% specify data: matrix of filenames and TR
%===========================================================================

clear tmp
for ses=1:nses
	dir{ses} = '../Data';		% directories can be generalised for >1 ses/sub
	tmp{ses} = spm_get('files',dir{ses},'swau*.img');
end
SPM.xY.P = cat(1,tmp{:});

SPM.xY.RT          = TR;


% Configure design matrix
%===========================================================================
SPM = spm_fmri_spm_ui(SPM);


% Estimate parameters
%===========================================================================
SPM = spm_spm(SPM);


% Add extra contrasts
%===========================================================================

ocon = length(SPM.xCon);

eb = eye(SPM.xBF.order);			% basis functions

clear cnam

cnam{1} = 'Faces vs Baseline';
cwgt{1} = [kron([1 0 0 0],eb) zeros(3,6)];		% zeros = movement params
ctyp{1} = 'F';

cnam{2} = 'Effects of lag';
cwgt{2} = [kron([0 1 0 0],eb) zeros(3,6)];	
ctyp{2} = 'F';

for c = 1:length(cnam)

	cw = [cwgt{c} zeros(size(cwgt{c},1),1)]';	% pad with zero for constant

	SPM.xCon(c+ocon)     = spm_FcUtil('Set',cnam{c},ctyp{c},'c',cw,SPM.xX.xKXs);
end


% and evaluate
%---------------------------------------------------------------------------
spm_contrasts(SPM);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now press RESULTS (and see README)!
% Parametric plots do not work for some reason?