Homec4science
Diffusion iCAPs public version 2 (master) Browse

iCAPs_v2/packages/spm12/config/spm_run_fmri_est.m
37b12c85c481master

Last Change
Raw File

spm_run_fmri_est.m
View Options

function out = spm_run_fmri_est(job)
% Estimate parameters of a specified model
% 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-2017 Wellcome Trust Centre for Neuroimaging
% $Id: spm_run_fmri_est.m 7057 2017-04-13 16:45:49Z guillaume $
%-Load SPM.mat file
%--------------------------------------------------------------------------
SPM = [];
load(job.spmmat{:});
out.spmmat = job.spmmat;
%-Move to the directory where the SPM.mat file is
%--------------------------------------------------------------------------
original_dir = pwd;
cd(fileparts(job.spmmat{:}));
%==========================================================================
% R E M L E S T I M A T I O N
%==========================================================================
if isfield(job.method,'Classical')
%-ReML estimation of the model
%----------------------------------------------------------------------
SPM = spm_spm(SPM);
%-Automatically set up contrasts for factorial designs
%----------------------------------------------------------------------
if isfield(SPM,'factor') && ~isempty(SPM.factor) && SPM.factor(1).levels > 1
Ic = [];
%-Generate contrasts
%------------------------------------------------------------------
cons = spm_design_contrasts(SPM);
%-Create F-contrasts
%------------------------------------------------------------------
for i=1:length(cons)
con = cons(i).c;
name = cons(i).name;
STAT = 'F';
[c,I] = spm_conman('ParseCon',con,SPM.xX.xKXs,STAT);
if all(I)
DxCon = spm_FcUtil('Set',name,STAT,'c',c,SPM.xX.xKXs);
if isempty(SPM.xCon)
SPM.xCon = DxCon;
else
SPM.xCon(end+1) = DxCon;
end
Ic = [Ic length(SPM.xCon)];
end
end
%-Create t-contrasts
%------------------------------------------------------------------
for i=1:length(cons)
% Create a t-contrast for each row of each F-contrast
% The resulting contrast image can be used in a 2nd-level analysis
Fcon = cons(i).c;
STAT = 'T';
for r=1:size(Fcon,1)
con = Fcon(r,:);
str = cons(i).name;
if strncmp('Interaction',str,11)
name = ['Positive ',str,'_',int2str(r)];
else
sp1 = find(isspace(str), 1);
name = ['Positive',str(sp1:end),'_',int2str(r)];
end
[c,I] = spm_conman('ParseCon',con,SPM.xX.xKXs,STAT);
if all(I)
DxCon = spm_FcUtil('Set',name,STAT,'c',c,SPM.xX.xKXs);
if isempty(SPM.xCon)
SPM.xCon = DxCon;
else
SPM.xCon(end+1) = DxCon;
end
Ic = [Ic length(SPM.xCon)];
end
end
end
%-Estimate constrasts
%------------------------------------------------------------------
if ~isempty(Ic)
spm('FnBanner','spm_contrasts.m');
SPM = spm_contrasts(SPM,Ic);
end
end
%-Residuals
%----------------------------------------------------------------------
if job.write_residuals
VRes = spm_write_residuals(SPM,NaN);
end
%-Computation results
%----------------------------------------------------------------------
%out.spmvar = SPM;
out.beta = spm_file({SPM.Vbeta(:).fname}','path',SPM.swd);
out.mask = {fullfile(SPM.swd,SPM.VM.fname)};
out.resms = {fullfile(SPM.swd,SPM.VResMS.fname)};
if job.write_residuals
out.res = spm_file({VRes(:).fname}','path',SPM.swd);
end
cd(original_dir);
fprintf('Done\n');
return
end
if job.write_residuals
warning('Save residuals option is only implemented for classical inference.');
job.write_residuals = false;
end
%==========================================================================
% B A Y E S I A N 2nd L E V E L E S T I M A T I O N
%==========================================================================
if isfield(job.method,'Bayesian2')
SPM = spm_spm_Bayes(SPM);
%out.spmvar = SPM;
cd(original_dir);
fprintf('Done\n');
return
end
%==========================================================================
% B A Y E S I A N 1st L E V E L E S T I M A T I O N
%==========================================================================
%-Analyse specific slices or whole volume
%--------------------------------------------------------------------------
switch char(fieldnames(job.method.Bayesian.space))
case 'volume'
SPM.PPM.space_type = 'volume';
SPM.PPM.block_type = lower(job.method.Bayesian.space.volume.block_type);
case 'slices'
SPM.PPM.space_type = 'slices';
SPM.PPM.AN_slices = job.method.Bayesian.space.slices.numbers;
SPM.PPM.block_type = lower(job.method.Bayesian.space.slices.block_type);
case 'clusters'
SPM.PPM.space_type = 'clusters';
SPM.PPM.clustermask = job.method.Bayesian.space.clusters.mask;
SPM.PPM.block_type = lower(job.method.Bayesian.space.clusters.block_type);
otherwise
SPM.PPM.space_type = 'volume';
SPM.PPM.block_type = 'slices';
end
%-Regression coefficient priors
%--------------------------------------------------------------------------
switch job.method.Bayesian.signal
case 'UGL'
SPM.PPM.priors.W = 'Spatial - UGL';
case 'GMRF'
SPM.PPM.priors.W = 'Spatial - GMRF';
case 'LORETA'
SPM.PPM.priors.W = 'Spatial - LORETA';
case 'WGL'
SPM.PPM.priors.W = 'Spatial - WGL';
case 'Global'
SPM.PPM.priors.W = 'Voxel - Shrinkage';
case 'Uninformative'
SPM.PPM.priors.W = 'Voxel - Uninformative';
otherwise
error('Unknown prior for W in Bayesian 1st level estimation.');
end
%-Number of AR coefficients
%--------------------------------------------------------------------------
SPM.PPM.AR_P = job.method.Bayesian.ARP;
%-AR coefficient priors
%--------------------------------------------------------------------------
if isfield(job.method.Bayesian.noise,'UGL')
SPM.PPM.priors.A = 'Spatial - UGL';
elseif isfield(job.method.Bayesian.noise,'GMRF')
SPM.PPM.priors.A = 'Spatial - GMRF';
elseif isfield(job.method.Bayesian.noise,'LORETA')
SPM.PPM.priors.A = 'Spatial - LORETA';
elseif isfield(job.method.Bayesian.noise,'tissue_type')
SPM.PPM.priors.A = 'Discrete';
SPM.PPM.priors.SY = char(job.method.Bayesian.noise.tissue_type);
elseif isfield(job.method.Bayesian.noise,'Robust')
SPM.PPM.priors.A = 'Robust';
SPM.PPM.AR_P = 0;
SPM.PPM.update_F = 1;
end
%-Define an empty contrast
%--------------------------------------------------------------------------
NullCon = spm_FcUtil('Set','','P','c',[],1);
NullCon.X0 = [];
NullCon.iX0 = [];
NullCon.X1o = [];
NullCon.eidf = 1;
SPM.xCon = [];
%-Set up contrasts for 2nd-level ANOVA
%--------------------------------------------------------------------------
if strcmp(job.method.Bayesian.anova.second,'Yes') && isfield(SPM,'factor')
cons = spm_design_contrasts(SPM);
for i=1:length(cons)
% Create a simple contrast for each row of each F-contrast
% The resulting contrast image can be used in a 2nd-level analysis
Fcon = cons(i).c;
for r=1:size(Fcon,1)
% Normalise contrast st. sum of positive elements is 1
% and sum of negative elements is 1
con = Fcon(r,:);
con = con / length(find(con==1));
% Change name
str = cons(i).name;
sp1 = find(str==' ', 1);
if strcmp(str(1:11),'Interaction')
name = ['Positive ',str,'_',int2str(r)];
else
name = ['Positive',str(sp1:end),'_',int2str(r)];
end
DxCon = NullCon;
DxCon.name = name;
DxCon.c = con';
if isempty(SPM.xCon)
SPM.xCon = DxCon;
else
SPM.xCon(end+1) = DxCon;
end
end
end
end
%-Set up user-specified simple contrasts
%--------------------------------------------------------------------------
K = size(SPM.xX.X,2);
for c = 1:length(job.method.Bayesian.gcon)
DxCon = NullCon;
DxCon.name = job.method.Bayesian.gcon(c).name;
convec = job.method.Bayesian.gcon(c).convec(:);
if length(convec) == K
DxCon.c = convec;
elseif length(convec) < K
fprintf('Zero padding.'); %-#
DxCon.c = [convec; zeros(K-length(convec),1)];
else
warning(['User-specified contrast nb %d has %d entries '...
'but there are %d regressors - ignored.'], c,length(convec),K);
DxCon = [];
end
if isempty(SPM.xCon)
SPM.xCon = DxCon;
else
SPM.xCon(end+1) = DxCon;
end
end
%-Compute log evidence maps
%--------------------------------------------------------------------------
if strcmp(job.method.Bayesian.LogEv,'Yes')
SPM.PPM.update_F = 1;
SPM.PPM.compute_det_D = 1;
end
%-1st level Bayesian ANOVA ?
%--------------------------------------------------------------------------
bayes_anova = 0;
if strcmp(job.method.Bayesian.anova.first,'Yes')
bayes_anova = 1;
SPM.PPM.update_F = 1; % Compute evidence for each model
SPM.PPM.compute_det_D = 1;
end
%-Variational Bayes estimation
%--------------------------------------------------------------------------
SPM = spm_spm_vb(SPM);
%-Bayesian ANOVA using model comparison
%--------------------------------------------------------------------------
if bayes_anova
% We don't want to estimate contrasts for each different model
SPM.xCon = [];
spm_vb_ppm_anova(SPM);
end
%out.spmvar = SPM;
cd(original_dir);
fprintf('Done\n')
c4science · Help