123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260 |
- function R = spm_searchlight(SPM,searchopt,fun,varargin)
- % Local mass-multivariate (c.f., searchlight) facility
- % FORMAT R = spm_searchlight(SPM,searchopt,fun,varargin)
- % SPM - structure with fields:
- % .xY.VY - filenames char array or spm_vol struct array of images
- % .VM - filename or spm_vol structure to a mask (binary) image
- % Mask image can have any orientation, voxel size or data type.
- % It is interpolated using nearest neighbour interpolation to
- % the voxel locations of the data.
- % If empty, all voxels are used.
- % searchopt - searchlight options using VOI structure (xY) from spm_ROI
- % .def - searchlight definition {['sphere'] 'box'}
- % .spec - searchlight parameters [sphere radius {mm}]
- % fun - function handle to a function that takes three input arguments:
- % a [n x v] matrix (nb images x nb voxels within searchlight)
- % a [3 x v] matrix of voxels location within searchlight {vox}
- % a list of parameters provided in varargin
- % and returns a vector value [1 x N]
- % varargin - list of parameters sent to fun
- %
- % R - a [N x 1] cell array with each output (fun nargout) reshaped
- % to a volume or directly a volume if N == 1
- % Values outside the mask are attributed NaN.
- %__________________________________________________________________________
- %
- % References:
- %
- % [1] Adaptive Analysis of fMRI Data. Friman O, Borga M, Lundberg P and
- % Knutsson H. (2003) NeuroImage 19(3):837-845.
- %
- % [2] Information-based functional brain mapping. Kriegeskorte N, Goebel R,
- % Bandettini P. (2006) PNAS 103: 3863-3868.
- %__________________________________________________________________________
- % Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging
- % Guillaume Flandin
- % $Id: spm_searchlight.m 4475 2011-09-09 17:53:14Z guillaume $
- spm('FnBanner',mfilename); %-#
- %-Get input images
- %--------------------------------------------------------------------------
- try
- VY = SPM.xY.VY;
- catch
- [VY, sts] = spm_select([1 Inf],'image','Select images');
- if ~sts, R = {}; return; end
- end
- if iscellstr(VY), VY = char(VY); end
- if ~isstruct(VY)
- VY = spm_vol(VY);
- end
- %-Check dimensions and orientations of all images
- %--------------------------------------------------------------------------
- spm_check_orientations(VY);
- %-Get mask image
- %--------------------------------------------------------------------------
- try
- VM = SPM.VM;
- catch
- [VM, sts] = spm_select([0 1],'image','Select mask');
- if ~sts, VM = []; end
- end
- if ~isstruct(VM) && ~isempty(VM)
- VM = spm_vol(VM);
- end
- %-Get space details
- %--------------------------------------------------------------------------
- N = numel(VY); %-number of images
- M = VY(1).mat; %-voxels to mm matrix
- iM = inv(M); %-mm to voxels matrix
- DIM = VY(1).dim; %-image dimensions
- NDIM = prod([DIM N]); %-overall dimension
- [x,y,z] = ndgrid(1:DIM(1),1:DIM(2),1:DIM(3));
- XYZ = [x(:)';y(:)';z(:)']; clear x y z %-voxel coordinates {vx}
- XYZmm = M(1:3,:)*[XYZ; ones(1,size(XYZ,2))];%-voxel coordinates {mm}
- XYZmm_cpy = XYZmm; %-copy without masking
- %-Strategy for reading data: memory [1] / block [2] / disk [3]
- %--------------------------------------------------------------------------
- fprintf('%-40s: ','Read data'); %-#
- try
- ds = NDIM*8; %-data size {bytes}
- MAXMEM = spm('memory'); %-max amnt usable memory
- if ds > MAXMEM
- YY = struct('y',[], 'i',[Inf Inf]);
- blk = floor(MAXMEM/(prod([DIM(1:2) N])*8));%-blk size {# slices}
- if blk == 0, error('revert to disk.'); end
- Ystrtg = 2;
- fprintf('%30s\n','...per block'); %-#
- else
- YY = spm_read_vols(VY);
- Ystrtg = 1;
- fprintf('%30s\n','...in memory'); %-#
- end
- catch
- Ystrtg = 3;
- fprintf('%30s\n','...from disk'); %-#
- end
- %-Search volume (from mask)
- %--------------------------------------------------------------------------
- fprintf('%-40s: ','Read mask'); %-#
- if ~isempty(VM)
- if any(DIM-VM.dim) || any(any(abs(VM.mat-M)>1e-4))
- MM = spm_get_data(VM,VM.mat\[XYZmm;ones(1,size(XYZmm,2))],false);
- else
- MM = spm_read_vols(VM);
- end
- MM = logical(MM);
- XYZmm = XYZmm(:,MM(:));
- XYZ = XYZ(:,MM(:));
- fprintf('%30s\n','...done'); %-#
- else
- MM = true(DIM);
- fprintf('%30s\n', '...none'); %-#
- end
- %-Searchlight options (clique definition)
- %--------------------------------------------------------------------------
- try, xY = searchopt; end
- xY.xyz = [NaN NaN NaN];
- xY.rej = {'cluster','mask'};
- xY = spm_ROI(xY);
- %-Evaluated function
- %--------------------------------------------------------------------------
- if ischar(fun)
- fun = str2func(fun);
- end
- if ~isa(fun, 'function_handle')
- error('''fun'' must be a function handle with two input parameters.');
- end
- %-Get local clique and perform searchlight over voxels
- %==========================================================================
- %-Build local clique
- %--------------------------------------------------------------------------
- fprintf('%-40s: ','Construct clique'); %-#
- c = round(DIM(:)/2);
- xY.xyz = M(1:3,:) * [c;1];
- [xY, clique] = spm_ROI(xY,XYZmm_cpy);
- clique = round(iM(1:3,:) * [clique;ones(1,size(clique,2))]);
- clique = bsxfun(@minus, clique, c);
- dc = (max(clique,[],2) - min(clique,[],2) + 1)';
- fprintf('%30s\n',sprintf('%d voxels - [%dx%dx%d]',size(clique,2),dc)); %-#
- if Ystrtg == 2 && blk < dc(3)
- fprintf('%-40s: %30s\n','Read data','revert to disk'); %-#
- Ystrtg = 3;
- end
- %-Initialise progress bar
- %--------------------------------------------------------------------------
- spm_figure('GetWin','Interactive');
- spm_progress_bar('Init',size(XYZ,2));
- Ibar = floor(linspace(1,size(XYZ,2), 100));
- fprintf('%-40s: %30s','Searchlight','...computing'); %-#
- SLR = [];
- %-Searchlight
- %--------------------------------------------------------------------------
- for i=1:size(XYZ,2)
- %-Local clique (handle image boundaries and mask)
- %----------------------------------------------------------------------
- xyz = bsxfun(@plus,XYZ(:,i),clique);
- xyz(:,any(bsxfun(@lt,xyz,[1 1 1]') | bsxfun(@gt,xyz,DIM'))) = [];
- idx = sub2ind(DIM,xyz(1,:),xyz(2,:),xyz(3,:));
- j = MM(idx);
- idx = idx(j);
- xyz = xyz(:,j);
-
- %-Read data
- %----------------------------------------------------------------------
- if Ystrtg == 3
- Y = spm_get_data(VY,xyz,false);
- elseif Ystrtg == 2
- if min(xyz(3,:)) < min(YY.i) || max(xyz(3,:)) > max(YY.i)
- YY.i = min(xyz(3,:)):min(xyz(3,:))+blk;
- YY.i(YY.i>DIM(3)) = [];
- YY.y = zeros(DIM(1),DIM(2),numel(YY.i),N);
- for v=1:N
- for p=1:numel(YY.i)
- YY.y(:,:,p,v) = ...
- spm_slice_vol(VY(v), ...
- spm_matrix([0 0 YY.i(p)]), ...
- VY(v).dim(1:2),0);
- end
- end
- end
- idx = idx - (min(YY.i)-1) * prod(DIM(1:2));
- k = prod(DIM(1:2))*numel(YY.i);
- idx = bsxfun(@plus, idx(:), 0:k:k*N-1);
- Y = YY.y(idx)';
- elseif Ystrtg == 1
- idx = bsxfun(@plus, idx(:), 0:prod(DIM):NDIM-1);
- Y = YY(idx)';
- end
-
- %-Call user-specified function
- %----------------------------------------------------------------------
- if isempty(SLR)
- t = fun(Y,xyz,varargin{:});
- SLR = zeros(size(XYZ,2),length(t));
- SLR(1,:) = t;
- else
- SLR(i,:) = fun(Y,xyz,varargin{:});
- end
-
- %-Update progress bar
- %----------------------------------------------------------------------
- if any(Ibar == i), spm_progress_bar('Set', i); end
-
- end
- %-Clear progress bar
- %--------------------------------------------------------------------------
- fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done'); %-#
- spm_progress_bar('Clear');
- %-Return computations - reshaped as volumes in a cell array
- %--------------------------------------------------------------------------
- R = cell(size(SLR,2),1);
- for i=1:size(SLR,2)
- MV = NaN(DIM);
- MV(sub2ind(DIM,XYZ(1,:),XYZ(2,:),XYZ(3,:))) = SLR(:,i);
- R{i} = MV;
- end
- %-Write images if required
- %--------------------------------------------------------------------------
- fprintf('%-40s: %30s','Output images','...writing'); %-#
- VO(1:size(SLR,2)) = deal(struct(...
- 'fname', [],...
- 'dim', DIM,...
- 'dt', [spm_type('float64') spm_platform('bigend')],...
- 'mat', M,...
- 'pinfo', [1 0 0]',...
- 'descrip', ['spm_searchlight: ' func2str(fun)]));
- for i=1:size(SLR,2)
- VO(i).fname = sprintf('%s_%04d%s', ...
- 'searchlight',i,spm_file_ext);
- VO(i).descrip = sprintf('%s (%04d)',VO(i).descrip,i);
- end
- for i=1:size(SLR,2)
- spm_write_vol(VO(i),R{i});
- end
- fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done'); %-#
- %-And exit
- %--------------------------------------------------------------------------
- if size(SLR,2) == 1, R = R{1}; end
- fprintf('%-40s: %30s\n','Completed',spm('time')); %-#
|