function [xVi, mask] = spm_est_non_sphericity(SPM) % Non-sphericity estimation using ReML % FORMAT [xVi, mask] = spm_est_non_sphericity(SPM) % % Required fields of SPM structure (see spm_spm): % SPM.xY.VY - nScan x 1 struct array of file handles % SPM.xX - structure containing design matrix information % SPM.xX.W - optional whitening/weighting matrix % SPM.xVi - structure describing intrinsic non-sphericity % SPM.xM - structure containing masking information % % Return xVi from SPM.xVi with extra fields: % xVi.V - estimated non-sphericity, trace(V) = rank(V) % xVi.h - hyperparameters xVi.V = xVi.h(1)*xVi.Vi{1} + ... % xVi.Cy - spatially whitened (used by ReML to estimate h) % % mask - logical array of voxels within analysis mask %__________________________________________________________________________ % % In a first pass, voxels over which non-sphericity will be estimated are % selected using an 'effects of interest' F-contrast (can be specified in % SPM.xVi.Fcontrast) and critical threshold taken from SPM defaults % stats..UFp. % The sample covariance matrix (xVi.Cy) is then estimated by pooling over % these voxels, assuming V is constant over them. % Finally, SPM will invoke ReML to estimate hyperparameters (xVi.h) of an % array of non-sphericity components (xVi.Vi), providing a high precise % estimate of the non-sphericity matrix (xVi.V). %__________________________________________________________________________ % Copyright (C) 1994-2016 Wellcome Trust Centre for Neuroimaging % Karl Friston & Guillaume Flandin % $Id: spm_est_non_sphericity.m 6913 2016-10-31 10:12:27Z guillaume $ SVNid = '$Rev: 6913 $'; %-Say hello %-------------------------------------------------------------------------- spm('FnBanner',mfilename,SVNid); %-Get data, design, mask and variance components details %-------------------------------------------------------------------------- VY = SPM.xY.VY; M = VY(1).mat; DIM = VY(1).dim; YNaNrep = spm_type(VY(1).dt(1),'nanrep'); xX = SPM.xX; nScan = size(xX.X,1); if ~isfield(xX,'W') xX.W = speye(nScan,nScan); end xM = SPM.xM; xVi = SPM.xVi; %-Compute Hsqr and F-threshold under i.i.d. %-------------------------------------------------------------------------- xX.xKXs = spm_sp('Set',spm_filter(xX.K,xX.W*xX.X)); xX.xKXs.X = full(xX.xKXs.X); xX.pKX = spm_sp('x-',xX.xKXs); if isfield(xVi,'Fcontrast') Fcname = 'User-specified contrast'; xCon = spm_FcUtil('Set',Fcname,'F','c',xVi.Fcontrast,xX.xKXs); else Fcname = 'effects of interest'; iX0 = [xX.iB xX.iG]; xCon = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs); end if ~isempty(xCon(1).c) X1o = spm_FcUtil('X1o', xCon(1),xX.xKXs); Hsqr = spm_FcUtil('Hsqr',xCon(1),xX.xKXs); trMV = spm_SpUtil('trMV',X1o); else % Force all voxels to enter non-sphericity trMV = 1; Hsqr = Inf; end trRV = spm_SpUtil('trRV',xX.xKXs); %-Threshold for voxels entering non-sphericity estimates %-------------------------------------------------------------------------- try modality = lower(spm_get_defaults('modality')); UFp = spm_get_defaults(['stats.' modality '.ufp']); catch UFp = 0.001; end xVi.UFp = UFp; UF = spm_invFcdf(1 - UFp,[trMV,trRV]); %========================================================================== %- P O O L E D V A R I A N C E E S T I M A T I O N %========================================================================== %-Get explicit mask(s) %========================================================================== mask = true(DIM); for i = 1:numel(xM.VM) if ~(isfield(SPM,'xVol') && isfield(SPM.xVol,'G')) %-Assume it fits entirely in memory C = spm_bsplinc(xM.VM(i), [0 0 0 0 0 0]'); v = true(DIM); [x1,x2] = ndgrid(1:DIM(1),1:DIM(2)); for x3 = 1:DIM(3) M2 = inv(M\xM.VM(i).mat); y1 = M2(1,1)*x1+M2(1,2)*x2+(M2(1,3)*x3+M2(1,4)); y2 = M2(2,1)*x1+M2(2,2)*x2+(M2(2,3)*x3+M2(2,4)); y3 = M2(3,1)*x1+M2(3,2)*x2+(M2(3,3)*x3+M2(3,4)); v(:,:,x3) = spm_bsplins(C, y1,y2,y3, [0 0 0 0 0 0]') > 0; end mask = mask & v; clear C v x1 x2 x3 M2 y1 y2 y3 else if spm_mesh_detect(xM.VM(i)) v = xM.VM(i).private.cdata() > 0; else v = spm_mesh_project(gifti(SPM.xVol.G), xM.VM(i)) > 0; end mask = mask & v(:); clear v end end Cy = 0; %- spatially whitened Cm = mask; %-mask of selected voxels %-Split data into 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,'Hyperparameter estimation','Chunks'); for i=1:nbchunks chunk = chunks(i):chunks(i+1)-1; %-Report progress %====================================================================== if i > 1, fprintf(repmat(sprintf('\b'),1,72)); end %-# fprintf('%-40s: %30s', sprintf('Chunk %3d/%-3d',i,nbchunks),... '...processing'); %-# %-Get data & construct analysis mask %---------------------------------------------------------------------- Y = zeros(nScan,numel(chunk)); cmask = mask(chunk); for j=1:nScan if ~any(cmask), break, end %-Break if empty mask Y(j,cmask) = spm_data_read(VY(j),chunk(cmask));%-Read chunk of data cmask(cmask) = Y(j,cmask) > xM.TH(j); %-Threshold (& NaN) mask if xM.I && ~YNaNrep && xM.TH(j) < 0 %-Use implicit mask cmask(cmask) = abs(Y(j,cmask)) > eps; end end cmask(cmask) = any(diff(Y(:,cmask),1)); %-Mask constant data mask(chunk) = cmask; Cm(chunk) = cmask; if ~any(cmask), continue, end Y = Y(:,cmask); %-Data within mask %-Remove filter confounds %---------------------------------------------------------------------- KWY = spm_filter(xX.K,xX.W*Y); %-Ordinary Least Squares estimation %---------------------------------------------------------------------- beta = xX.pKX*KWY; %-Parameter estimates if any(cmask) res = spm_sp('r',xX.xKXs,KWY); %-Residuals else res = zeros(nScan,0); end ResSS = sum(res.^2); %-Residual SSQ clear res %-F-threshold & accumulate spatially whitened Y*Y' %---------------------------------------------------------------------- j = sum((Hsqr*beta).^2,1)/trMV > UF*ResSS/trRV; Cm(chunk(cmask)) = j; q = nnz(j); if q q = spdiags(sqrt(trRV./ResSS(j)'),0,q,q); Y = Y(:,j)*q; Cy = Cy + Y*Y'; end %-Report progress %====================================================================== fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...done'); %-# spm_progress_bar('Set',i); end fprintf('\n'); %-# spm_progress_bar('Clear'); s = nnz(Cm); %-Number of selected voxels if ~s error('Please check your data: There are no significant voxels.'); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %-Estimate sample covariance matrix Y*Y' from voxels in Cm % Cy = 0; % nbchunks = ceil(s / chunksize); % chunks = min(cumsum([1 repmat(chunksize,1,nbchunks)]),s+1); % vx = find(Cm); % for i=1:nbchunks % chunk = chunks(i):chunks(i+1)-1; % n = numel(chunk); % Y = zeros(nScan,n); % for j=1:nScan % Y(j,:) = spm_data_read(VY(j),vx(chunk)); % end % % store ResSS above to prevent recomputing it here: % %-------------------------------------------------- % KY = spm_filter(xX.K,Y); % ResSS = spm_sp('r',xX.xKXs,KY); % ResSS = sum(ResSS.^2); % %-------------------------------------------------- % q = spdiags(sqrt(trRV./ResSS'),0,n,n); % Y = Y*q; % Cy = Cy + Y*Y'; % end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Cy = Cy / s; %-Sample covariance matrix %========================================================================== %- R e M L E S T I M A T I O N %========================================================================== %-ReML estimate of residual correlations through hyperparameters (h) %-------------------------------------------------------------------------- str = sprintf('Temporal non-sphericity (%d voxels)',s); fprintf('%-40s: %30s\n',str,'...ReML estimation'); %-# %-ReML for separable designs and covariance components %-------------------------------------------------------------------------- if isstruct(xX.K) m = length(xVi.Vi); h = zeros(m,1); V = sparse(nScan,nScan); for i = 1:length(xX.K) % extract blocks from bases %------------------------------------------------------------------ q = xX.K(i).row; p = []; Qp = {}; for j = 1:m if nnz(xVi.Vi{j}(q,q)) Qp{end + 1} = xVi.Vi{j}(q,q); p = [p j]; end end % design space for ReML (with confounds in filter) %------------------------------------------------------------------ Xp = xX.X(q,:); try Xp = [Xp xX.K(i).X0]; end % ReML %------------------------------------------------------------------ fprintf('%-30s\n',sprintf(' ReML Block %i',i)); [Vp,hp] = spm_reml(Cy(q,q),Xp,Qp); V(q,q) = V(q,q) + Vp; h(p) = hp; end else [V,h] = spm_reml(Cy,xX.X,xVi.Vi); end %-Normalize non-sphericity and save hyperparameters %-------------------------------------------------------------------------- V = V*nScan/trace(V); xVi.h = h; xVi.V = V; % Save non-sphericity xVi.V xVi.Cy = Cy; % spatially whitened