Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/Utilities/spm12/spm_write_residuals.m
8a186cd82567master

Last Change
Raw File

spm_write_residuals.m
View Options

function VRes = spm_write_residuals(SPM,Ic)
% Write residual images
% FORMAT Vres = spm_write_residuals(SPM,Ic)
% SPM - structure containing generic analysis details
% Ic - contrast index used to adjust data (0: no adjustment)
% (NaN: adjust for everything)
%
% VRes - struct array of residual image handles
%__________________________________________________________________________
% Copyright (C) 2012-2013 Wellcome Trust Centre for Neuroimaging
% Guillaume Flandin
% $Id: spm_write_residuals.m 5575 2013-07-08 16:38:51Z guillaume $
%-Get SPM.mat
%--------------------------------------------------------------------------
if ~nargin || isempty(SPM)
[SPM,sts] = spm_select(1,'^SPM\.mat$','Select SPM.mat');
if ~sts, VRes = ''; return; end
end
if ~isstruct(SPM)
swd = spm_file(SPM,'fpath');
try
load(fullfile(swd,'SPM.mat'));
SPM.swd = swd;
catch
error(['Cannot read ' fullfile(swd,'SPM.mat')]);
end
end
try, SPM.swd; catch, SPM.swd = pwd; end
cwd = pwd; cd(SPM.swd);
%-Get contrast used to adjust data
%--------------------------------------------------------------------------
if nargin<2 || isempty(Ic)
q(1) = 0;
Con = {'<don''t adjust>'};
q(2) = NaN;
Con{2} = '<adjust for everything>';
for i = 1:length(SPM.xCon)
if strcmp(SPM.xCon(i).STAT,'F')
q(end + 1) = i;
Con{end + 1} = SPM.xCon(i).name;
end
end
i = spm_input('adjust data for (select contrast)','!+1','m',Con);
Ic = q(i);
end
%-Compute and write residuals
%==========================================================================
spm('Pointer','Watch')
M = SPM.xY.VY(1).mat;
DIM = SPM.xY.VY(1).dim(1:3);
[nScan, nBeta] = size(SPM.xX.X);
if spm_mesh_detect(SPM.xY.VY)
file_ext = '.gii';
else
file_ext = spm_file_ext;
end
%-Initialise residual images
%--------------------------------------------------------------------------
VRes(1:nScan) = deal(struct(...
'fname', [],...
'dim', DIM,...
'dt', [spm_type('float64') spm_platform('bigend')],...
'mat', M,...
'pinfo', [1 0 0]',...
'descrip', 'Residuals'));
for i = 1:nScan
VRes(i).fname = [sprintf('Res_%04d', i) file_ext];
VRes(i).descrip = sprintf('Residuals (%04d)', i);
end
VRes = spm_data_hdr_write(VRes);
%-Loop over chunks
%--------------------------------------------------------------------------
chunksize = floor(spm_get_defaults('stats.maxmem') / 8 / nScan);
nbchunks = ceil(prod(DIM) / chunksize);
chunks = min(cumsum([1 repmat(chunksize,1,nbchunks)]),prod(DIM)+1);
spm_progress_bar('Init',nbchunks,'Writing residuals','Chunks');
for i=1:nbchunks
chunk = chunks(i):chunks(i+1)-1;
%-Get mask
%----------------------------------------------------------------------
m = spm_data_read(SPM.VM,chunk) > 0;
m = m(:)';
%-Get raw data, whiten and filter
%----------------------------------------------------------------------
y = zeros(nScan,numel(chunk));
for j=1:nScan
y(j,:) = spm_data_read(SPM.xY.VY(j),chunk);
end
y(:,~m) = [];
y = spm_filter(SPM.xX.K,SPM.xX.W*y);
if Ic ~= 0
%-Parameter estimates: beta = xX.pKX*xX.K*y
%------------------------------------------------------------------
beta = zeros(nBeta,numel(chunk));
for j=1:nBeta
beta(j,:) = spm_data_read(SPM.Vbeta(j),chunk);
end
beta(:,~m) = [];
%-Subtract Y0 = XO*beta, Y = Yc + Y0 + e
%------------------------------------------------------------------
if ~isnan(Ic)
y = y - spm_FcUtil('Y0',SPM.xCon(Ic),SPM.xX.xKXs,beta);
else
y = y - SPM.xX.xKXs.X * beta;
end
end
%-Write residuals
%----------------------------------------------------------------------
yy = NaN(numel(chunk),1);
for j=1:nScan
yy(m) = y(j,:);
VRes(j) = spm_data_write(VRes(j), yy, chunk);
end
spm_progress_bar('Set',i)
end
cd(cwd);
spm_progress_bar('Clear')
spm('Pointer','Arrow')
c4science · Help