1 function mb_5_slice_timing(P, sliceorder, refslice, timing, imethod, prefix)
   2 % function mb_5_slice_timing(P, sliceorder, refslice, timing, imethod, prefix)
   3 % INPUT:
   4 % 	P		nimages x ?	Matrix with filenames
   5 %                       can also be a cell array of the above (multiple subj).
   6 %	sliceorder	slice acquisition order, a vector of integers, each
   7 %                       integer referring the slice number in the image file
   8 %			(1=first), and the order of integers representing their
   9 %			temporal acquisition order
  10 %	refslice	slice for time 0
  11 %	timing		additional information for sequence timing
  12 %			timing(1) = time between slices
  13 %			timing(2) = time between last slices and next
  14 %			volume
  15 %       imethod          (see help interp1) one of ['sinc']
  16 %                       'sinc'    - using custom function here
  17 %                       'nearest' - nearest neighbor interpolation
  18 %                       'linear'  - linear interpolation
  19 %                       'spline'  - cubic spline interpolation
  20 %                       'cubic'   - cubic interpolation
  21 %       prefix          prefix for new file names ['a']
  22 %		      
  23 % 	If no input is specified the function serves as a GUI			
  24 %
  25 % OUTPUT:
  26 %	None
  27 %
  28 %   Slice-time corrected files are prepended with an 'a'
  29 %
  30 %   Note: The sliceorder arg that specifies slice acquisition order is
  31 %   a vector of N numbers, where N is the number of slices per volume.
  32 %   Each number refers to the position of a slice within the image file.
  33 %   The order of numbers within the vector is the temporal order in which
  34 %   those slices were acquired.
  35 %
  36 %   To check the order of slices within an image file, use the SPM Display
  37 %   option and move the crosshairs to a voxel co-ordinate of z=1.  This
  38 %   corresponds to a point in the first slice of the volume.
  39 %
  40 %   The function corrects differences in slice acquisition times.
  41 %   This routine is intended to correct for the staggered order of
  42 %   slice acquisition that is used during echoplanar scanning. The
  43 %   correction is necessary to make the data on each slice correspond
  44 %   to the same point in time. Without correction, the data on one
  45 %   slice will represent a point in time as far removed as 1/2 the TR
  46 %   from an adjacent slice (in the case of an interleaved sequence).
  47 %
  48 %   This routine "shifts" a signal in time to provide an output
  49 %   vector that represents the same (continuous) signal sampled
  50 %   starting either later or earlier. This is accomplished by a simple
  51 %   shift of the phase of the sines that make up the signal.
  52 %
  53 %   Recall that a Fourier transform allows for a representation of any
  54 %   signal as the linear combination of sinusoids of different
  55 %   frequencies and phases. Effectively, we will add a constant
  56 %   to the phase of every frequency, shifting the data in time.
  57 %
  58 %    Shifter - This is the filter by which the signal will be convolved
  59 %    to introduce the phase shift. It is constructed explicitly in
  60 %    the Fourier domain. In the time domain, it may be described as
  61 %    an impulse (delta function) that has been shifted in time the
  62 %    amount described by TimeShift.
  63 %
  64 %   The correction works by lagging (shifting forward) the time-series
  65 %     data on each slice using sinc-interpolation. This results in each
  66 %     time series having the values that would have been obtained had
  67 %     the slice been acquired at the same time as the reference slice.
  68 %
  69 %   To make this clear, consider a neural event (and ensuing hemodynamic
  70 %     response) that occurs simultaneously on two adjacent slices. Values
  71 %     from slice "A" are acquired starting at time zero, simultaneous to
  72 %     the neural event, while values from slice "B" are acquired one
  73 %     second later. Without corection, the "B" values will describe a
  74 %     hemodynamic response that will appear to have began one second
  75 %     EARLIER on the "B" slice than on slice "A". To correct for this,
  76 %     the "B" values need to be shifted towards the Right, i.e., towards
  77 %     the last value.
  78 %
  79 %   This correction assumes that the data are band-limited (i.e. there
  80 %     is no meaningful information present in the data at a frequency
  81 %     higher than that of the Nyquist). This assumption is support by
  82 %     the study of Josephs et al (1997, NeuroImage) that obtained
  83 %     event-related data at an effective TR of 166 msecs. No physio-
  84 %     logical signal change was present at frequencies higher than our
  85 %     typical Nyquist (0.25 HZ).
  86 %
  87 % Written by Darren Gitelman at Northwestern U., 1998
  88 %
  89 % Based (in large part) on ACQCORRECT.PRO from Geoff Aguirre and
  90 % Eric Zarahn at U. Penn.
  91 %
  92 % v1.0	07/04/98	DRG
  93 % v1.1  07/09/98	DRG	fixed code to reflect 1-based indices
  94 %				of matlab vs. 0-based of pvwave
  95 %
  96 % Modified by R Henson, C Buechel and J Ashburner, FIL, to
  97 % handle different reference slices and memory mapping.
  98 %
  99 % Modified by M Erb, at U. Tuebingen, 1999, to ask for non-continuous
 100 % slice timing and number of sessions.
 101 %
 102 % Modified by R Henson for more general slice order and SPM2
 103 %_______________________________________________________________________
 104 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
 105 
 106 %
 107 % $Id: spm_slice_timing.m 671 2006-11-02 12:08:04Z john $
 108 
 109 
 110 SPMid = spm('FnBanner',mfilename,'$Rev: 671 $');
 111 [Finter,Fgraph,CmdLine] = spm('FnUIsetup','Slice timing');
 112 spm_help('!ContextHelp',mfilename);
 113 
 114 if nargin < 1,
 115         % get number of subjects
 116         nsubjects = spm_input('number of subjects/sessions',1, 'e', 1);
 117         if nsubjects < 1,
 118                 spm_figure('Clear','Interactive');
 119                 return;
 120         end
 121 
 122 	for i = 1:nsubjects,
 123 		% Choose the images
 124 		PP = [];
 125 		PP = spm_select(Inf,'image',...
 126 			['Select images to acquisition correct for subject ' num2str(i)]);
 127 		P{i} = PP;
 128 	end;
 129 end;
 130 
 131 if iscell(P),
 132 	nsubjects = length(P);
 133 else,
 134 	nsubjects = 1;
 135 	P = {P};
 136 end;
 137 
 138 Vin 	= spm_vol(P{1}(1,:));
 139 nslices	= Vin(1).dim(3);
 140 
 141 % Modified (simplified) by R. Henson	03/06/25
 142 if nargin < 2,
 143   	sliceorder=[];
 144 	while length(sliceorder)~=nslices | max(sliceorder)>nslices | ...
 145 		min(sliceorder)<1 | any(diff(sort(sliceorder))~=1),
 146 		sliceorder = spm_input(...
 147 		   'Acquisition order? (1=first slice in image)','!+0','e');
 148 	end;
 149 end;
 150 % End of modification by R. Henson
 151 
 152 if nargin < 3,
 153 	% Choose reference slice (in Analyze format, slice 1 = bottom)
 154 	% Note: no checking that 1 < refslice < no.slices (default = middle slice)
 155 	refslice = spm_input('Reference slice? (1=first slice in image)','!+0','e',floor(Vin(1).dim(3)/2));
 156 end;
 157 
 158 if nargin < 4,
 159 	TR = spm_input('Interscan interval (TR) {secs}','!+1','e',3);
 160 	TA = spm_input('Acquisition Time (TA) {secs}','!+1','e',TR-TR/nslices);
 161 	while TA > TR | TA <= 0,
 162 		TA = spm_input('Acquisition Time (TA) {secs}','!+0','e',TA);
 163 	end;
 164 	timing(2) = TR - TA;
 165 	timing(1) = TA / (nslices -1);
 166 	factor = timing(1)/TR;
 167 else,
 168 	TR 	= (nslices-1)*timing(1)+timing(2);
 169 	fprintf('Your TR is %1.1f\n',TR);
 170 	factor = timing(1)/TR;
 171 end;
 172 if nargin < 5
 173   imethod = 'sinc';
 174 end
 175 if nargin < 6
 176   prefix = 'a';
 177 end
 178 
 179 spm('Pointer','Watch')
 180 
 181 for subj = 1:nsubjects
 182 	task = ['Slice timing: working on session ' num2str(subj)];
 183 	spm('FigName',task,Finter,CmdLine);
 184 	PP      = P{subj};
 185 	Vin 	= spm_vol(PP);
 186 	nimgo	= numel(Vin);
 187 	if strcmp(imethod, 'sinc')
 188 	  nimg	= 2^(floor(log2(nimgo))+1);
 189 	else
 190 	  nimg = nimgo;
 191 	end 
 192 	nslices_t= Vin(1).dim(3);
 193 	if nslices_t ~= nslices,
 194 		fprintf('Number of slices differ! %d %\n', nimg);
 195 	else
 196 
 197 		% create new header files
 198 		Vout 	= Vin;
 199 		for k=1:nimgo,
 200 			[pth,nm,xt,vr] = fileparts(deblank(Vin(k).fname));
 201 			Vout(k).fname  = fullfile(pth,[prefix nm xt vr]);
 202 			if isfield(Vout(k),'descrip'),
 203 				desc = [Vout(k).descrip ' '];
 204 			else,
 205 				desc = '';
 206 			end;
 207 			Vout(k).descrip = [desc 'acq-fix ref-slice ' int2str(refslice)];
 208 		end;
 209 		Vout = spm_create_vol(Vout);
 210 
 211 		% Set up large matrix for holding image info
 212 		% Organization is time by voxels
 213 		slices = zeros([Vout(1).dim(1:2) nimgo]);
 214 		stack  = zeros([nimg Vout(1).dim(1)]);
 215 
 216 		spm_progress_bar('Init',nslices,'Correcting acquisition delay','planes complete');
 217 
 218 		% For loop to read data slice by slice do correction and write out
 219 		% In analzye format, the first slice in is the first one in the volume.
 220 
 221 		rslice = find(sliceorder==refslice);
 222 		for k = 1:nslices,
 223 
 224 		  % Set up time acquired within slice order
 225 		  shiftamount  = (find(sliceorder==k) - rslice) * factor;
 226 		  
 227 		  % Read in slice data
 228 		  B  = spm_matrix([0 0 k]);
 229 		  for m=1:nimgo,
 230 		    slices(:,:,m) = spm_slice_vol(Vin(m),B,Vin(1).dim(1:2),1);
 231 		  end;
 232 		  
 233 		  if strcmp(imethod, 'sinc')
 234 		    
 235 		    % set up shifting variables
 236 		    len     = size(stack,1);
 237 		    phi     = zeros(1,len);
 238 		    
 239 		    % Check if signal is odd or even -- impacts how Phi is reflected
 240 		    %  across the Nyquist frequency. Opposite to use in pvwave.
 241 		    OffSet  = 0;
 242 		    if rem(len,2) ~= 0, OffSet = 1; end;
 243 		    
 244 		    % Phi represents a range of phases up to the Nyquist frequency
 245 		    % Shifted phi 1 to right.
 246 		    for f = 1:len/2,
 247 		      phi(f+1) = -1*shiftamount*2*pi/(len/f);
 248 		    end;
 249 		    
 250 		    % Mirror phi about the center
 251 		    % 1 is added on both sides to reflect Matlab's 1 based indices
 252 		    % Offset is opposite to program in pvwave again because indices are 1 based
 253 		    phi(len/2+1+1-OffSet:len) = -fliplr(phi(1+1:len/2+OffSet));
 254 		    
 255 		    % Transform phi to the frequency domain and take the complex transpose
 256 		    shifter = [cos(phi) + sin(phi)*sqrt(-1)].';
 257 		    shifter = shifter(:,ones(size(stack,2),1)); % Tony's trick
 258 		    
 259 		  else % i.e. not sinc
 260 		    
 261 		    % interpolation vector - see help interp1
 262 		    X1 = (1:nimgo)-shiftamount;
 263 		    imeth = ['*' imethod];
 264 		  end
 265 		  % Loop over columns
 266 		  for i=1:Vout(1).dim(2),
 267 		    
 268 		    % Extract columns from slices
 269 		    stack(1:nimgo,:) = reshape(slices(:,i,:),[Vout(1).dim(1) nimgo])';
 270 		    
 271 		    if strcmp(imethod, 'sinc')
 272 		      
 273 		      % fill in continous function to avoid edge effects
 274 		      for g=1:size(stack,2),
 275 			stack(nimgo+1:end,g) = linspace(stack(nimgo,g),...
 276 							stack(1,g),nimg-nimgo)';
 277 		      end;
 278 		      
 279 		      % shift the columns
 280 		      stack = real(ifft(fft(stack,[],1).* ...
 281 					shifter,[],1));
 282 		    else
 283 		      % interpolate within columns
 284 		      stack = interp1(stack, X1, imeth);
 285 		    end
 286 		    
 287 		    % Re-insert shifted columns
 288 		    slices(:,i,:) = reshape(stack(1:nimgo,:)',[Vout(1).dim(1) 1 nimgo]);
 289 		  end;
 290 		  
 291 		  % write out the slice for all volumes
 292 		  for p = 1:nimgo,
 293 		    Vout(p) = spm_write_plane(Vout(p),slices(:,:,p),k);
 294 		  end;
 295 		  spm_progress_bar('Set',k);
 296 		end;
 297 		spm_progress_bar('Clear');
 298 	end
 299 end
 300 
 301 spm('FigName','Slice timing: done',Finter,CmdLine);
 302 spm('Pointer');
 303 return;