function mb_5_slice_timing(P, sliceorder, refslice, timing, imethod, prefix)
% function mb_5_slice_timing(P, sliceorder, refslice, timing, imethod, prefix)
% INPUT:
% P nimages x ? Matrix with filenames
% can also be a cell array of the above (multiple subj).
% sliceorder slice acquisition order, a vector of integers, each
% integer referring the slice number in the image file
% (1=first), and the order of integers representing their
% temporal acquisition order
% refslice slice for time 0
% timing additional information for sequence timing
% timing(1) = time between slices
% timing(2) = time between last slices and next
% volume
% imethod (see help interp1) one of ['sinc']
% 'sinc' - using custom function here
% 'nearest' - nearest neighbor interpolation
% 'linear' - linear interpolation
% 'spline' - cubic spline interpolation
% 'cubic' - cubic interpolation
% prefix prefix for new file names ['a']
%
% If no input is specified the function serves as a GUI
%
% OUTPUT:
% None
%
% Slice-time corrected files are prepended with an 'a'
%
% Note: The sliceorder arg that specifies slice acquisition order is
% a vector of N numbers, where N is the number of slices per volume.
% Each number refers to the position of a slice within the image file.
% The order of numbers within the vector is the temporal order in which
% those slices were acquired.
%
% To check the order of slices within an image file, use the SPM Display
% option and move the crosshairs to a voxel co-ordinate of z=1. This
% corresponds to a point in the first slice of the volume.
%
% The function corrects differences in slice acquisition times.
% This routine is intended to correct for the staggered order of
% slice acquisition that is used during echoplanar scanning. The
% correction is necessary to make the data on each slice correspond
% to the same point in time. Without correction, the data on one
% slice will represent a point in time as far removed as 1/2 the TR
% from an adjacent slice (in the case of an interleaved sequence).
%
% This routine "shifts" a signal in time to provide an output
% vector that represents the same (continuous) signal sampled
% starting either later or earlier. This is accomplished by a simple
% shift of the phase of the sines that make up the signal.
%
% Recall that a Fourier transform allows for a representation of any
% signal as the linear combination of sinusoids of different
% frequencies and phases. Effectively, we will add a constant
% to the phase of every frequency, shifting the data in time.
%
% Shifter - This is the filter by which the signal will be convolved
% to introduce the phase shift. It is constructed explicitly in
% the Fourier domain. In the time domain, it may be described as
% an impulse (delta function) that has been shifted in time the
% amount described by TimeShift.
%
% The correction works by lagging (shifting forward) the time-series
% data on each slice using sinc-interpolation. This results in each
% time series having the values that would have been obtained had
% the slice been acquired at the same time as the reference slice.
%
% To make this clear, consider a neural event (and ensuing hemodynamic
% response) that occurs simultaneously on two adjacent slices. Values
% from slice "A" are acquired starting at time zero, simultaneous to
% the neural event, while values from slice "B" are acquired one
% second later. Without corection, the "B" values will describe a
% hemodynamic response that will appear to have began one second
% EARLIER on the "B" slice than on slice "A". To correct for this,
% the "B" values need to be shifted towards the Right, i.e., towards
% the last value.
%
% This correction assumes that the data are band-limited (i.e. there
% is no meaningful information present in the data at a frequency
% higher than that of the Nyquist). This assumption is support by
% the study of Josephs et al (1997, NeuroImage) that obtained
% event-related data at an effective TR of 166 msecs. No physio-
% logical signal change was present at frequencies higher than our
% typical Nyquist (0.25 HZ).
%
% Written by Darren Gitelman at Northwestern U., 1998
%
% Based (in large part) on ACQCORRECT.PRO from Geoff Aguirre and
% Eric Zarahn at U. Penn.
%
% v1.0 07/04/98 DRG
% v1.1 07/09/98 DRG fixed code to reflect 1-based indices
% of matlab vs. 0-based of pvwave
%
% Modified by R Henson, C Buechel and J Ashburner, FIL, to
% handle different reference slices and memory mapping.
%
% Modified by M Erb, at U. Tuebingen, 1999, to ask for non-continuous
% slice timing and number of sessions.
%
% Modified by R Henson for more general slice order and SPM2
%_______________________________________________________________________
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
%
% $Id: spm_slice_timing.m 671 2006-11-02 12:08:04Z john $
SPMid = spm('FnBanner',mfilename,'$Rev: 671 $');
[Finter,Fgraph,CmdLine] = spm('FnUIsetup','Slice timing');
spm_help('!ContextHelp',mfilename);
if nargin < 1,
% get number of subjects
nsubjects = spm_input('number of subjects/sessions',1, 'e', 1);
if nsubjects < 1,
spm_figure('Clear','Interactive');
return;
end
for i = 1:nsubjects,
% Choose the images
PP = [];
PP = spm_select(Inf,'image',...
['Select images to acquisition correct for subject ' num2str(i)]);
P{i} = PP;
end;
end;
if iscell(P),
nsubjects = length(P);
else,
nsubjects = 1;
P = {P};
end;
Vin = spm_vol(P{1}(1,:));
nslices = Vin(1).dim(3);
% Modified (simplified) by R. Henson 03/06/25
if nargin < 2,
sliceorder=[];
while length(sliceorder)~=nslices | max(sliceorder)>nslices | ...
min(sliceorder)<1 | any(diff(sort(sliceorder))~=1),
sliceorder = spm_input(...
'Acquisition order? (1=first slice in image)','!+0','e');
end;
end;
% End of modification by R. Henson
if nargin < 3,
% Choose reference slice (in Analyze format, slice 1 = bottom)
% Note: no checking that 1 < refslice < no.slices (default = middle slice)
refslice = spm_input('Reference slice? (1=first slice in image)','!+0','e',floor(Vin(1).dim(3)/2));
end;
if nargin < 4,
TR = spm_input('Interscan interval (TR) {secs}','!+1','e',3);
TA = spm_input('Acquisition Time (TA) {secs}','!+1','e',TR-TR/nslices);
while TA > TR | TA <= 0,
TA = spm_input('Acquisition Time (TA) {secs}','!+0','e',TA);
end;
timing(2) = TR - TA;
timing(1) = TA / (nslices -1);
factor = timing(1)/TR;
else,
TR = (nslices-1)*timing(1)+timing(2);
fprintf('Your TR is %1.1f\n',TR);
factor = timing(1)/TR;
end;
if nargin < 5
imethod = 'sinc';
end
if nargin < 6
prefix = 'a';
end
spm('Pointer','Watch')
for subj = 1:nsubjects
task = ['Slice timing: working on session ' num2str(subj)];
spm('FigName',task,Finter,CmdLine);
PP = P{subj};
Vin = spm_vol(PP);
nimgo = numel(Vin);
if strcmp(imethod, 'sinc')
nimg = 2^(floor(log2(nimgo))+1);
else
nimg = nimgo;
end
nslices_t= Vin(1).dim(3);
if nslices_t ~= nslices,
fprintf('Number of slices differ! %d %\n', nimg);
else
% create new header files
Vout = Vin;
for k=1:nimgo,
[pth,nm,xt,vr] = fileparts(deblank(Vin(k).fname));
Vout(k).fname = fullfile(pth,[prefix nm xt vr]);
if isfield(Vout(k),'descrip'),
desc = [Vout(k).descrip ' '];
else,
desc = '';
end;
Vout(k).descrip = [desc 'acq-fix ref-slice ' int2str(refslice)];
end;
Vout = spm_create_vol(Vout);
% Set up large matrix for holding image info
% Organization is time by voxels
slices = zeros([Vout(1).dim(1:2) nimgo]);
stack = zeros([nimg Vout(1).dim(1)]);
spm_progress_bar('Init',nslices,'Correcting acquisition delay','planes complete');
% For loop to read data slice by slice do correction and write out
% In analzye format, the first slice in is the first one in the volume.
rslice = find(sliceorder==refslice);
for k = 1:nslices,
% Set up time acquired within slice order
shiftamount = (find(sliceorder==k) - rslice) * factor;
% Read in slice data
B = spm_matrix([0 0 k]);
for m=1:nimgo,
slices(:,:,m) = spm_slice_vol(Vin(m),B,Vin(1).dim(1:2),1);
end;
if strcmp(imethod, 'sinc')
% set up shifting variables
len = size(stack,1);
phi = zeros(1,len);
% Check if signal is odd or even -- impacts how Phi is reflected
% across the Nyquist frequency. Opposite to use in pvwave.
OffSet = 0;
if rem(len,2) ~= 0, OffSet = 1; end;
% Phi represents a range of phases up to the Nyquist frequency
% Shifted phi 1 to right.
for f = 1:len/2,
phi(f+1) = -1*shiftamount*2*pi/(len/f);
end;
% Mirror phi about the center
% 1 is added on both sides to reflect Matlab's 1 based indices
% Offset is opposite to program in pvwave again because indices are 1 based
phi(len/2+1+1-OffSet:len) = -fliplr(phi(1+1:len/2+OffSet));
% Transform phi to the frequency domain and take the complex transpose
shifter = [cos(phi) + sin(phi)*sqrt(-1)].';
shifter = shifter(:,ones(size(stack,2),1)); % Tony's trick
else % i.e. not sinc
% interpolation vector - see help interp1
X1 = (1:nimgo)-shiftamount;
imeth = ['*' imethod];
end
% Loop over columns
for i=1:Vout(1).dim(2),
% Extract columns from slices
stack(1:nimgo,:) = reshape(slices(:,i,:),[Vout(1).dim(1) nimgo])';
if strcmp(imethod, 'sinc')
% fill in continous function to avoid edge effects
for g=1:size(stack,2),
stack(nimgo+1:end,g) = linspace(stack(nimgo,g),...
stack(1,g),nimg-nimgo)';
end;
% shift the columns
stack = real(ifft(fft(stack,[],1).* ...
shifter,[],1));
else
% interpolate within columns
stack = interp1(stack, X1, imeth);
end
% Re-insert shifted columns
slices(:,i,:) = reshape(stack(1:nimgo,:)',[Vout(1).dim(1) 1 nimgo]);
end;
% write out the slice for all volumes
for p = 1:nimgo,
Vout(p) = spm_write_plane(Vout(p),slices(:,:,p),k);
end;
spm_progress_bar('Set',k);
end;
spm_progress_bar('Clear');
end
end
spm('FigName','Slice timing: done',Finter,CmdLine);
spm('Pointer');
return;