123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119 |
- function Vo = spm_resss(Vi,Vo,R,flags)
- % Create residual sum of squares image (ResSS)
- % FORMAT Vo = spm_resss(Vi,Vo,R,flags)
- % Vi - vector of mapped image volumes to work on (from spm_vol)
- % Vo - handle structure for mapped output image volume
- % R - residual forming matrix
- % flags - 'm' for implicit zero masking
- % Vo (output) - handle structure of output image volume after modifications
- % for writing
- %
- % Note that spm_create_vol needs to be called external to this function -
- % the header is not created.
- %__________________________________________________________________________
- %
- % Residuals are computed as R*Y, where Y is the data vector read from
- % images mapped as Vi. The residual sum of squares image (mapped as Vo)
- % is written.
- %
- %--------------------------------------------------------------------------
- %
- % For a simple linear model Y = X*B * E, with design matrix X,
- % (unknown) parameter vector(s) B, and data matrix Y, the least squares
- % estimates of B are given by b = inv(X'*X)*X'*Y. If X is rank
- % deficient, then the Moore-Penrose pseudoinverse may be used to obtain
- % the least squares parameter estimates with the lowest L2 norm: b =
- % pinv(X)*Y.
- %
- % The fitted values are then y = X*b = X*inv(X'*X)*X'*Y, (or
- % y=X*pinv(X)*Y). Since the fitted values y are usually known as
- % "y-hat", X*inv(X'*X)*X' is known as the "hat matrix" for this model,
- % denoted H.
- %
- % The residuals for this fit (estimates of E) are e = Y - y.
- % Substituting from the above, e = (I-H)*Y, where I is the identity
- % matrix (see eye). (I-H) is called the residual forming matrix,
- % denoted R.
- %
- % Geometrically, R is a projection matrix, projecting the data into the
- % subspace orthogonal to the design space.
- %
- % ----------------
- %
- % For temporally smoothed fMRI models with convolution matrix K, R is a
- % little more complicated:
- % K*Y = K*X * B + K*E
- % KY = KX * B + KE
- % ...a little working shows that hat matrix is H = KX*inv(KX'*KX)*KX'
- % (or KX*pinv(KX)), where KX=K*X. The smoothed residuals KE (=K*E) are
- % then given from the temporally smoothed data KY (=K*Y) by y=H*KY.
- % Thus the residualising matrix for the temporally smoothed residuals
- % from the temporally smoothed data is then (I-H).
- %
- % Usually the image time series is not temporally smoothed, in which
- % case the hat and residualising matrices must incorporate the temporal
- % smoothing. The hat matrix for the *raw* (unsmoothed) time series Y is
- % H*K, and the corresponding residualising matrix is R=(K-H*K).
- % In full, that's
- % R = (K - KX*inv(KX'*KX)*KX'*K)
- % or R = (K - KX*pinv(KX)*K) when using a pseudoinverse
- %
- %--------------------------------------------------------------------------
- %
- % This function can also be used when the b's are images. The residuals
- % are then e = Y - X*b, so let Vi refer to the vector of images and
- % parameter estimates ([Y;b]), and then R is ([eye(n),-X]), where n is
- % the number of Y images.
- %
- %--------------------------------------------------------------------------
- %
- % Don't forget to either apply any image scaling (grand mean or
- % proportional scaling global normalisation) to the image scalefactors,
- % or to combine the global scaling factors in the residual forming
- % matrix.
- %__________________________________________________________________________
- % Copyright (C) 1999-2012 Wellcome Trust Centre for Neuroimaging
- % Andrew Holmes & John Ashburner
- % $Id: spm_resss.m 5097 2012-12-06 16:08:16Z guillaume $
- %-Argument checks
- %--------------------------------------------------------------------------
- if nargin<4, flags=''; end, if isempty(flags), flags='-'; end
- mask = any(flags=='m');
- if nargin<3, error('insufficient arguments'); end
- ni = size(R,2); %-ni = #images
- if ni~=numel(Vi), error('incompatible dimensions'); end
- if Vo.dt(1)~=spm_type('float32')
- error('only float output images supported');
- end
- %-Image dimension, orientation and voxel size checks
- %--------------------------------------------------------------------------
- spm_check_orientations([Vi Vo]);
- %==========================================================================
- % - C O M P U T A T I O N
- %==========================================================================
- Y = zeros([Vo.dim(1:2),ni]); %-PlaneStack data
- im = false(ni,1);
- for j=1:ni
- im(j)=~spm_type(Vi(j).dt(1),'NaNrep'); %-Images without NaNrep
- end
- %-Loop over planes computing ResSS
- for p=1:Vo.dim(3)
- M = spm_matrix([0 0 p]); %-Sampling matrix
- %-Read plane data
- for j=1:ni, Y(:,:,j) = spm_slice_vol(Vi(j),M,Vi(j).dim(1:2),0); end
- %-Apply implicit zero mask for image types without a NaNrep
- if mask, Y(Y(:,:,im)==0)=NaN; end
- e = R*reshape(Y,prod(Vi(1).dim(1:2)),ni)'; %-residuals as DataMtx
- ss = reshape(sum(e.^2,1),Vi(1).dim(1:2)); %-ResSS plane
- Vo = spm_write_plane(Vo,ss,p); %-Write plane
- end
|