Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/packages/spm12/spm_contrasts.m
31022b3a13e6master

Last Change
Raw File

spm_contrasts.m
View Options

function SPM = spm_contrasts(SPM,Ic)
% Compute and store contrast parameters and inference SPM{.}
% FORMAT SPM = spm_contrasts(SPM,Ic)
%
% SPM - SPM data structure
% Ic - indices of xCon to compute
%
% This function fills in SPM.xCon and writes con_????, ess_???? and
% spm?_???? images.
%__________________________________________________________________________
% Copyright (C) 2002-2017 Wellcome Trust Centre for Neuroimaging
% Karl Friston, Will Penny & Guillaume Flandin
% $Id: spm_contrasts.m 7029 2017-02-24 15:39:07Z guillaume $
% Temporary copy of the SPM variable, to avoid saving it in SPM.mat unless
% it has changed (faster, read-only access)
%--------------------------------------------------------------------------
tmpSPM = SPM;
%-Change to results directory
%--------------------------------------------------------------------------
try, cd(SPM.swd); end
%-Get contrast definitions (if available)
%--------------------------------------------------------------------------
try
xCon = SPM.xCon;
catch
xCon = [];
end
%-Set all contrasts by default
%--------------------------------------------------------------------------
if nargin < 2
Ic = 1:length(xCon);
end
Ic(Ic == 0) = [];
%-Map parameter and hyperarameter files
%--------------------------------------------------------------------------
if ~isempty(xCon) && xCon(Ic(1)).STAT == 'P'
%-Conditional estimators
%----------------------------------------------------------------------
Vbeta = SPM.VCbeta;
else
%-OLS estimators and error variance estimate
%----------------------------------------------------------------------
Vbeta = SPM.Vbeta;
VHp = SPM.VResMS;
end
if spm_mesh_detect(Vbeta)
file_ext = '.gii';
g = SPM.xY.VY(1).private;
metadata = g.private.metadata;
name = {metadata.name};
if any(ismember(name,'SurfaceID'))
metadata = metadata(ismember(name,'SurfaceID'));
metadata = {metadata.name, metadata.value};
elseif isfield(g,'faces') && ~isempty(g.faces)
metadata = {'SurfaceID', SPM.xY.VY(1).fname};
else
metadata = {};
end
else
file_ext = spm_file_ext;
metadata = {};
end
%-Compute & store contrast parameters, contrast/ESS images, & SPM images
%==========================================================================
spm('Pointer','Watch')
XYZ = SPM.xVol.XYZ;
iXYZ = cumprod([1,SPM.xVol.DIM(1:2)'])*XYZ - sum(cumprod(SPM.xVol.DIM(1:2)'));
for i = 1:length(Ic)
%-Canonicalise contrast structure with required fields
%----------------------------------------------------------------------
ic = Ic(i);
if isempty(xCon(ic).eidf)
X1o = spm_FcUtil('X1o',xCon(ic),SPM.xX.xKXs);
[trMV,trMVMV] = spm_SpUtil('trMV',X1o,SPM.xX.V);
xCon(ic).eidf = trMV^2/trMVMV;
end
%-Write contrast/ESS images?
%======================================================================
if isempty(xCon(ic).Vcon)
switch xCon(ic).STAT
case {'T','P'}
if strcmp(xCon(ic).STAT,'P') && strcmp(SPM.PPM.xCon(ic).PSTAT,'F')
% Bayes Factor for compound contrast
%------------------------------------------------------
disp('Bayes factor for compound contrast');
fprintf('\t%-32s: %30s',sprintf('LogBF image %2d',ic),...
'...computing'); %-#
if isfield(SPM.PPM,'VB')
% First level Bayes
xCon = spm_vb_logbf(SPM,XYZ,xCon,ic);
else
% Second level Bayes
xCon = spm_bayes2_logbf(SPM,XYZ,xCon,ic);
end
else
%-Implement contrast as linear combination of beta images
%------------------------------------------------------
fprintf('\t%-32s: %30s',sprintf('contrast image %2d',ic),...
'...computing'); %-#
%-Prepare handle for contrast image
%------------------------------------------------------
xCon(ic).Vcon = struct(...
'fname', [sprintf('con_%04d',ic) file_ext],...
'dim', SPM.xVol.DIM',...
'dt', [spm_type('float32'), spm_platform('bigend')],...
'mat', SPM.xVol.M,...
'pinfo', [1,0,0]',...
'descrip',sprintf('Contrast %d: %s',ic,xCon(ic).name),...
metadata{:});
xCon(ic).Vcon = spm_data_hdr_write(xCon(ic).Vcon);
%-Compute contrast
%------------------------------------------------------
Q = find(abs(xCon(ic).c) > 0);
V = Vbeta(Q);
cB = zeros(1,size(XYZ,2));
for j=1:numel(V)
cB = cB + xCon(ic).c(Q(j)) * spm_data_read(V(j),'xyz',XYZ);
end
%-Write contrast image
%------------------------------------------------------
tmp = NaN(SPM.xVol.DIM');
tmp(iXYZ) = cB;
xCon(ic).Vcon = spm_data_write(xCon(ic).Vcon,tmp);
clear tmp cB
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),sprintf(...
'...written %s',spm_file(xCon(ic).Vcon.fname,'filename')))%-#
end
case 'F' %-Implement ESS as sum of squared weighted beta images
%----------------------------------------------------------
fprintf('\t%-32s: %30s',sprintf('ESS image %2d',ic),...
'...computing'); %-#
%-Prepare handle for ESS image
%----------------------------------------------------------
xCon(ic).Vcon = struct(...
'fname', [sprintf('ess_%04d',ic) file_ext],...
'dim', SPM.xVol.DIM',...
'dt', [spm_type('float32'), spm_platform('bigend')],...
'mat', SPM.xVol.M,...
'pinfo', [1,0,0]',...
'descrip',sprintf('ESS contrast %d: %s',ic,xCon(ic).name),...
metadata{:});
xCon(ic).Vcon = spm_data_hdr_write(xCon(ic).Vcon);
%-Compute ESS
%----------------------------------------------------------
% Residual (in parameter space) forming matrix
h = spm_FcUtil('Hsqr',xCon(ic),SPM.xX.xKXs);
ss = zeros(numel(Vbeta),size(XYZ,2));
for j=1:numel(Vbeta)
ss(j,:) = spm_data_read(Vbeta(j),'xyz',XYZ);
end
ss = sum((h*ss).^2,1);
%-Write ESS image
%----------------------------------------------------------
tmp = NaN(SPM.xVol.DIM');
tmp(iXYZ) = ss;
xCon(ic).Vcon = spm_data_write(xCon(ic).Vcon,tmp);
clear tmp ss
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),sprintf(...
'...written %s',spm_file(xCon(ic).Vcon.fname,'filename')))%-#
otherwise
%----------------------------------------------------------
error(['unknown STAT "',xCon(ic).STAT,'"'])
end % (switch(xCon...)
end % (if isempty(xCon(ic)...)
%-Write inference SPM/PPM
%======================================================================
if isempty(xCon(ic).Vspm) || xCon(ic).STAT == 'P'
% (always update PPM as size threshold, gamma, may have changed)
fprintf('\t%-32s: %30s',sprintf('spm{%s} image %2d',xCon(ic).STAT,ic),...
'...computing'); %-#
switch(xCon(ic).STAT)
case 'T' %-Compute SPM{t} image
%----------------------------------------------------------
cB = spm_data_read(xCon(ic).Vcon,'xyz',XYZ);
l = spm_data_read(VHp,'xyz',XYZ); % get hyperparamters
Vc = xCon(ic).c'*SPM.xX.Bcov*xCon(ic).c;
SE = sqrt(l*Vc); % and standard error
Z = cB./SE;
str = sprintf('[%.1f]',SPM.xX.erdf);
case 'P' %-Compute PPM{P} image
%----------------------------------------------------------
if all(strcmp({SPM.PPM.xCon(ic).PSTAT},'T'))
% Simple contrast - Gaussian distributed
c = xCon(ic).c;
cB = spm_data_read(xCon(ic).Vcon,'xyz',XYZ);
if isfield(SPM.PPM,'VB');
% If posterior sd image for that contrast does
% not already exist, then compute it
try
SPM.PPM.Vcon_sd(ic);
catch
SPM = spm_vb_contrasts(SPM,XYZ,xCon,ic);
end
% Read in posterior sd image for contrast
Vsd = spm_data_read(SPM.PPM.Vcon_sd(ic),'xyz',XYZ);
VcB = Vsd.^2;
else
VcB = c'*SPM.PPM.Cby*c;
for j = 1:length(SPM.PPM.l)
% hyperparameter and Taylor approximation
%----------------------------------------------
l = spm_data_read(SPM.VHp(j),'xyz',XYZ);
VcB = VcB + (c'*SPM.PPM.dC{j}*c)*(l - SPM.PPM.l(j));
end
end
% posterior probability cB > g
%------------------------------------------------------
Gamma = xCon(ic).eidf;
Z = 1 - spm_Ncdf(Gamma,cB,VcB);
% Convert probability to Log Odds Ratio
Z = log( Z ./ (1 - Z+eps) );
str = sprintf('[%.2f]',Gamma);
%xCon(ic).name = [xCon(ic).name ' ' str];
else
% Compound contrast - Log Bayes Factor
fprintf('\t\t%-75s\n','Log Bayes Factor for compound contrast');
fprintf('\t%-32s: %29s\n',' ',' ');
Z = spm_data_read(xCon(ic).Vcon,'xyz',XYZ);
str = sprintf('[%1.2f]',xCon(ic).eidf);
end
case 'F' %-Compute SPM{F} image
%----------------------------------------------------------
MVM = spm_data_read(xCon(ic).Vcon,'xyz',XYZ)/trMV;
RVR = spm_data_read(VHp,'xyz',XYZ);
Z = MVM./RVR;
str = sprintf('[%.1f,%.1f]',xCon(ic).eidf,SPM.xX.erdf);
otherwise
%----------------------------------------------------------
error(['unknown STAT "',xCon(ic).STAT,'"']);
end % (switch(xCon(ic)...)
%-Write SPM - statistic image
%------------------------------------------------------------------
xCon(ic).Vspm = struct(...
'fname', [sprintf('spm%s_%04d',xCon(ic).STAT,ic) file_ext],...
'dim', SPM.xVol.DIM',...
'dt', [spm_type('float32'), spm_platform('bigend')],...
'mat', SPM.xVol.M,...
'pinfo', [1,0,0]',...
'descrip',sprintf('SPM{%s_%s} - contrast %d: %s',...
xCon(ic).STAT,str,ic,xCon(ic).name),...
metadata{:});
xCon(ic).Vspm = spm_data_hdr_write(xCon(ic).Vspm);
tmp = zeros(SPM.xVol.DIM');
tmp(iXYZ) = Z;
xCon(ic).Vspm = spm_data_write(xCon(ic).Vspm,tmp);
clear tmp Z
cmd = sprintf(['[hReg,xSPM,SPM] = spm_results_ui(''Setup'',',...
'struct(''swd'',''%s'',''Ic'',%d));',...
'TabDat = spm_list(''List'',xSPM,hReg);'],pwd,ic);
img = spm_file(spm_file(xCon(ic).Vspm.fname,'filename'),'link',cmd);
n = 30; if length(img)>n, n = length(img)+n-13; end
fprintf('%s%*s\n',repmat(sprintf('\b'),1,30),n,sprintf(...
'...written %s',img)); %-#
end % (if isempty(xCon(ic)...)
end % (for i = 1:length(Ic))
spm('Pointer','Arrow')
% place xCon back in SPM
%--------------------------------------------------------------------------
SPM.xCon = xCon;
% Check if SPM has changed. Save only if it has.
%--------------------------------------------------------------------------
if spm_check_version('matlab','8.0') >= 0, my_isequaln = @isequaln;
else my_isequaln = @isequalwithequalnans; end
if ~my_isequaln(tmpSPM,SPM)
fprintf('\t%-32s: %30s','Saving SPM.mat','...writing'); %-#
save('SPM.mat', 'SPM', spm_get_defaults('mat.format'));
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...SPM.mat saved') %-#
end
c4science · Help