Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/packages/spm12/config/spm_run_fmri_spec.m
31022b3a13e6master

Last Change
Raw File

spm_run_fmri_spec.m
View Options

function out = spm_run_fmri_spec(job)
% Setting up the general linear model for fMRI time-series
% SPM job execution function
% takes a harvested job data structure and call SPM functions to perform
% computations on the data.
% Input:
% job - harvested job data structure (see matlabbatch help)
% Output:
% out - computation results, usually a struct variable.
%__________________________________________________________________________
% Copyright (C) 2005-2015 Wellcome Trust Centre for Neuroimaging
% $Id: spm_run_fmri_spec.m 6562 2015-09-25 14:38:17Z guillaume $
%-Check presence of previous analysis
%==========================================================================
%-Output directory
%--------------------------------------------------------------------------
cwd = pwd;
d = spm_file(job.dir{1},'cpath');
if ~exist(d,'dir')
sts = mkdir(d);
if ~sts, error('Error creating output directory "%s".',d); end
end
cd(d);
%-Ask about overwriting files from previous analyses
%--------------------------------------------------------------------------
if exist(fullfile(pwd,'SPM.mat'),'file')
str = {'Current directory contains existing SPM file:',...
'Continuing will overwrite existing file!'};
if spm_input(str,1,'bd','stop|continue',[1,0],1,mfilename);
fprintf('%-40s: %30s\n\n',...
'Abort... (existing SPM file)',spm('time'));
out = []; return
end
end
%-Delete old analysis files
%--------------------------------------------------------------------------
files = {'^mask\..{3}$','^ResMS\..{3}$','^RPV\..{3}$',...
'^beta_.{4}\..{3}$','^con_.{4}\..{3}$','^ResI_.{4}\..{3}$',...
'^ess_.{4}\..{3}$', '^spm\w{1}_.{4}\..{3}$'};
for i=1:length(files)
j = spm_select('List',pwd,files{i});
for k=1:size(j,1)
spm_unlink(deblank(j(k,:)));
end
end
%-Timing parameters & Basis functions
%==========================================================================
%-Repeat time
%--------------------------------------------------------------------------
SPM.xY.RT = job.timing.RT;
%-Basis function parameters
%--------------------------------------------------------------------------
SPM.xBF.UNITS = job.timing.units;
SPM.xBF.T = job.timing.fmri_t;
SPM.xBF.T0 = job.timing.fmri_t0;
%-Basis functions
%--------------------------------------------------------------------------
bf = char(fieldnames(job.bases));
if strcmp(bf,'none')
SPM.xBF.name = 'NONE';
elseif strcmp(bf,'hrf')
if all(job.bases.hrf.derivs == [0 0])
SPM.xBF.name = 'hrf';
elseif all(job.bases.hrf.derivs == [1 0])
SPM.xBF.name = 'hrf (with time derivative)';
elseif all(job.bases.hrf.derivs == [1 1])
SPM.xBF.name = 'hrf (with time and dispersion derivatives)';
else
error('Unknown HRF derivative choices.');
end
else
switch bf
case 'fourier'
SPM.xBF.name = 'Fourier set';
case 'fourier_han'
SPM.xBF.name = 'Fourier set (Hanning)';
case 'gamma'
SPM.xBF.name = 'Gamma functions';
case 'fir'
SPM.xBF.name = 'Finite Impulse Response';
otherwise
error('Unknown basis functions.');
end
SPM.xBF.length = job.bases.(bf).length;
SPM.xBF.order = job.bases.(bf).order;
end
%-Model interactions (Volterra)
%--------------------------------------------------------------------------
SPM.xBF.Volterra = job.volt;
%-Data & Design
%==========================================================================
design_only = ~isfield(job,'mask');
if ~design_only
SPM.xY.P = [];
end
for i = 1:numel(job.sess)
sess = job.sess(i);
%-Image filenames
%----------------------------------------------------------------------
if design_only
SPM.nscan(i) = sess.nscan;
else
SPM.nscan(i) = numel(sess.scans);
if SPM.nscan(i) == 1
sess.scans = cellstr(spm_select('Expand',sess.scans{1}));
SPM.nscan(i) = numel(sess.scans);
end
SPM.xY.P = strvcat(SPM.xY.P,sess.scans{:});
end
if SPM.nscan(i) == 1
error('Not enough scans in session %d.',i);
end
%-Multiple conditions (structure from a MAT-file)
%----------------------------------------------------------------------
if ~isempty(sess.multi{1})
%-Load MAT-file
%------------------------------------------------------------------
try
multicond = load(sess.multi{1});
catch
error('Cannot load %s',sess.multi{1});
end
%-Check structure content
%------------------------------------------------------------------
if ~all(isfield(multicond, {'names','onsets','durations'})) || ...
~iscell(multicond.names) || ...
~iscell(multicond.onsets) || ...
~iscell(multicond.durations) || ...
~isequal(numel(multicond.names), numel(multicond.onsets), ...
numel(multicond.durations))
error(['Multiple conditions MAT-file ''%s'' is invalid:\n',...
'File must contain names, onsets, and durations '...
'cell arrays of equal length.\n'],sess.multi{1});
end
for j=1:numel(multicond.onsets)
%-Mutiple Conditions: names, onsets and durations
%--------------------------------------------------------------
cond.name = multicond.names{j};
if isempty(cond.name)
error('MultiCond file: sess %d cond %d has no name.',i,j);
end
cond.onset = multicond.onsets{j};
if isempty(cond.onset)
error('MultiCond file: sess %d cond %d has no onset.',i,j);
end
cond.duration = multicond.durations{j};
if isempty(cond.onset)
error('MultiCond file: sess %d cond %d has no duration.',i,j);
end
%-Mutiple Conditions: Time Modulation
%--------------------------------------------------------------
if ~isfield(multicond,'tmod');
cond.tmod = 0;
else
try
cond.tmod = multicond.tmod{j};
catch
error('Error specifying time modulation.');
end
end
%-Mutiple Conditions: Parametric Modulation
%--------------------------------------------------------------
cond.pmod = [];
if isfield(multicond,'pmod')
try
%-Check if a PM is defined for that condition
if (j <= numel(multicond.pmod)) && ...
~isempty(multicond.pmod(j).name)
for ii = 1:numel(multicond.pmod(j).name)
cond.pmod(ii).name = multicond.pmod(j).name{ii};
cond.pmod(ii).param = multicond.pmod(j).param{ii};
cond.pmod(ii).poly = multicond.pmod(j).poly{ii};
end
end
catch
warning('Error specifying parametric modulation.');
rethrow(lasterror);
end
end
%-Mutiple Conditions: Orthogonalisation of Modulations
%--------------------------------------------------------------
if isfield(multicond,'orth') && (j <= numel(multicond.orth))
cond.orth = multicond.orth{j};
else
cond.orth = true;
end
%-Append to singly-specified conditions
%--------------------------------------------------------------
sess.cond(end+1) = cond;
end
end
%-Conditions
%----------------------------------------------------------------------
U = [];
for j = 1:numel(sess.cond)
%-Name, Onsets, Durations
%------------------------------------------------------------------
cond = sess.cond(j);
U(j).name = {cond.name};
U(j).ons = cond.onset(:);
U(j).dur = cond.duration(:);
U(j).orth = cond.orth;
if isempty(U(j).orth), U(j).orth = true; end
if length(U(j).dur) == 1
U(j).dur = repmat(U(j).dur,size(U(j).ons));
elseif numel(U(j).dur) ~= numel(U(j).ons)
error('Mismatch between number of onset and number of durations.');
end
%-Modulations
%------------------------------------------------------------------
P = [];
q1 = 0;
%-Time Modulation
switch job.timing.units
case 'secs'
sf = 1 / 60;
case 'scans'
sf = job.timing.RT / 60;
otherwise
error('Unknown unit "%s".',job.timing.units);
end
if cond.tmod > 0
P(1).name = 'time';
P(1).P = U(j).ons * sf;
P(1).h = cond.tmod;
q1 = 1;
end
%-Parametric Modulations
if ~isempty(cond.pmod)
for q = 1:numel(cond.pmod)
q1 = q1 + 1;
P(q1).name = cond.pmod(q).name;
P(q1).P = cond.pmod(q).param(:);
P(q1).h = cond.pmod(q).poly;
end
end
%-None
if isempty(P)
P.name = 'none';
P.h = 0;
end
U(j).P = P;
end
SPM.Sess(i).U = U;
%-User-specified regressors
%----------------------------------------------------------------------
C = [];
Cname = cell(1,numel(sess.regress));
for q = 1:numel(sess.regress)
Cname{q} = sess.regress(q).name;
if numel(sess.regress(q).val(:)) ~= SPM.nscan(i)
error('Length of regressor is not commensurate with data points.');
end
C = [C, sess.regress(q).val(:)];
end
%-Multiple regressors (from a TXT/MAT-file)
%----------------------------------------------------------------------
if ~isempty(sess.multi_reg{1})
for q=1:numel(sess.multi_reg)
tmp = load(sess.multi_reg{q});
names = {};
if isstruct(tmp) % .mat
if isfield(tmp,'R')
R = tmp.R;
if isfield(tmp,'names')
names = tmp.names;
end
elseif isfield(tmp,'xY');
R = tmp.xY.u;
names = {tmp.xY.name};
elseif isfield(tmp,'PPI')
R = [tmp.PPI.ppi tmp.PPI.Y tmp.PPI.P];
names = {...
['PPI Interaction: ' tmp.PPI.name],...
['Main Effect: ' tmp.PPI.xY.name ' BOLD'],...
'Main Effect: "psych" condition'};
else
error(['Variable ''R'' not found in multiple ' ...
'regressor file ''%s''.'], sess.multi_reg{q});
end
elseif isnumeric(tmp) % .txt
R = tmp;
end
if size(R,1) ~= SPM.nscan(i)
error('Length of regressor is not commensurate with data points.');
end
if ~isempty(names) && numel(names) ~= size(R,2)
warning('Mismatch between number of regressors and their names.');
names = {};
end
C = [C, R];
if numel(sess.multi_reg) == 1
nb_mult = '';
else
nb_mult = sprintf('_%d',q);
end
for j=1:size(R,2)
if isempty(names)
Cname{end+1} = sprintf('R%d%s',j,nb_mult);
else
Cname{end+1} = sprintf('%s%s',names{j},nb_mult);
end
end
end
end
SPM.Sess(i).C.C = C;
SPM.Sess(i).C.name = Cname;
end
%-Factorial design
%--------------------------------------------------------------------------
if ~isempty(job.fact)
for i=1:numel(job.fact)
SPM.factor(i).name = job.fact(i).name;
SPM.factor(i).levels = job.fact(i).levels;
end
if prod([SPM.factor.levels]) ~= numel(SPM.Sess(1).U)
error('Factors do not match conditions');
end
else
SPM.factor = [];
end
%-Globals
%--------------------------------------------------------------------------
SPM.xGX.iGXcalc = job.global;
%-Masking threshold
%--------------------------------------------------------------------------
SPM.xM.gMT = job.mthresh;
%-High Pass filter
%--------------------------------------------------------------------------
for i = 1:numel(job.sess),
SPM.xX.K(i).HParam = job.sess(i).hpf;
end
%-Autocorrelation
%--------------------------------------------------------------------------
SPM.xVi.form = job.cvi;
%-Setting up the GLM
%==========================================================================
SPM = spm_fmri_spm_ui(SPM);
%-Explicit mask
%--------------------------------------------------------------------------
if ~design_only
if ~isempty(job.mask{1})
SPM.xM.VM = spm_data_hdr_read(job.mask{1});
SPM.xM.xs.Masking = [SPM.xM.xs.Masking, '+explicit mask'];
end
end
%-Save SPM.mat
%--------------------------------------------------------------------------
fprintf('%-40s: ','Saving SPM configuration') %-#
save('SPM.mat','SPM', spm_get_defaults('mat.format'));
fprintf('%30s\n','...SPM.mat saved') %-#
fprintf('%-40s: %30s\n','Completed',spm('time')) %-#
out.spmmat{1} = fullfile(pwd, 'SPM.mat');
%-Change back directory
%--------------------------------------------------------------------------
cd(cwd);
c4science · Help