Dear List:
I am forwarding some comments by Rik Henson regarding the slice
timing issue. Rik just got back from holiday and he may have to go
back again (certainly if I bring up the slice timing issue again...
just kidding, :) ).
Although I probably don't have any intellectual right to comment on
what Rik says, that's never stopped me before. I agree with his
comments. Much of the confusion (at least for me) has stemmed from a
misunderstanding of the canonical interleaved sequence included with
spm99. In any case I agree that everything works as specified if you
use that sequence. I believe the code I sent that Rik has included
below provides a simpler approach to slice specification (simpler is
good). It is based on time. That is no matter what sequence is
selected the user would enter the desired time bin, and it would
figure out the correct slice. So entering 1 means slice 10 in a 10
slice descending sequence and slice 1 in an ascending sequence. If
you use the code please let me know if there is a problem (please
Lord let there not be a problem). Be careful of the carriage returns
when extracting from email. I'll post a clean copy on our ftp site
ftp.cnadc.northwestern.edu /pub/spm in a couple days.
As always, hope this doesn't hurt, and thanks Rik.
Darren
-------------------- begin Rik Henson's email ----------------
Darren/Andreas -
Apologies for not replying sooner - I have been away on holiday.
The concept behind slice order specification is actually very simple,
but there appears to be much room for confusion.
Let me try to summarize the development of spm_slice_timing...
When Christian and I wrote the front-end for the slice-timing
code you provided (adapted from Geoff and Eric's code), we made
the following simple assumptions:
1. The temporal order of slices would be coded by the
left-to-right order of slices in a vector
2. The slices would be referenced by the Analyze convention,
where 1=bottom slice.
If you follow these rules when entering the slice order (as user-
specified), there is no problem. Thus a six-slice descending sequence
in space (top slice acquired first) would be coded as:
6 5 4 3 2 1
I repeat: no errors will have been made if you followed these rules,
regardless what sequence you use.
The potential for confusion perhaps arose when we introduced two further
options. The first was the the addition of various default menu options
(descending, ascending, interleaved) that avoided the need for the user
to type in slice order by hand (for the most common sequences). The
second
was the default choice of a reference slice (particularly in relation
to interleaved sequences)....
1. Default options: Interleaving
================================
The ascending/descending options refer to acquisition order in space,
and are equivalent to typing 1:N or N:-1:1 respectively as user-
specified (where N=number of slices).
However, there are many possible interleaved sequences one can use:
Odd-even slices, top-to-middle, bottom-to-middle. Rather than trying to
offer all these, we chose one particular interleaving: top-to-middle.
For
a six-slice sequence, this would give:
6 3 5 2 4 1
Note for an odd-number of slices, the "odd" slice is placed first
(ie for a five-slice sequence: 3 5 2 4 1, rather than the alternative:
5 3 4 2 1).
Your scanner may however use a different sequence. In this case, you
should type the sliceorder by hand in the user-specified option.
Thus it is possible that you might have pressed interleaved thinking
it was appropriate for the particular interleaving used on your scanner
when it was not. Our apologies for not making this clearer in the
menu options.
2. Default reference slice
==========================
As detailed in many previous emails to this list, the default reference
slice (to which all other slices are "synchronised") offered by SPM is
the middle slice in space (in Analyze format).
Remember: if the temporal interpolation were perfect, the choice of
reference slice would not matter. However, the interpolation will alias
frequencies above the Nyquist limit, which may introduce appreciable
noise if your TR is long (eg ~>2s).
The amount of noise introduced will tend to increase the more the
timeseries is shifted (ie the futher away, in time, a slice is from
the reference slice), up to TR/2.
Because we use sequential acquisition (ascending/descending) here, our
rationale for offering the middle slice in space was that, because
the middle slice in space happens to coincide with the middle slice
in time for sequential acquistions, the maximum interpolation error
would then be "pushed" to the top and bottom slices, which are usually
at the extremes of the field-of-view and of least interest (eg the top
and bottom of the brain, which usually have the least grey-matter
voxels).
If however you are interested in a particular region within your
field-of-view (other than the middle), you might want to change
the reference slice to coincide with that region.
With interleaved sequences however, the choice of reference slice
is far from clear (on the above rationale). For example, with a
top-to-middle sequence, eg:
10 5 9 4 8 3 7 2 6 1
The middle slice(s) in space (5/6) are not the middle in time (8/3).
This is also true for odd-even sequences, eg:
10 8 6 4 2 9 7 5 3 1
Without an a priori "slice-of-interest", there seems no principled
reason for choosing the reference slice with interleaved sequences.
In this case, one might as well chose the first slice in time (10
in the above examples), to avoid the need to change the default
fMRI_T0 in the SPM stats stage (see point 3 below). Thus you must
think before accepting the default reference slice (which was
chosen for sequential acquisitions) when using interleaved acquisitions.
Other points to note:
3. Reference slice and reference time-bin
=========================================
As discussed in many previous emails, the value of fMRI_T0 in the
Defaults-Statistics-fMRI menu should correspond to your choice of
reference slice. fMRI_T0 tells SPM which of the fMRI_T time-bins
(per TR) should be used as the value of the covariate for that scan.
Namely, if you have interpolated to slice r of N (where r now refers
to time, not Analyse format), then:
fMRI_T0 = fMRI_T * (r/N)
Note that the default value is 1 (ie first slice acquired). This makes
sense for epoch designs (eg boxcars), where slice-timing is less of a
problem.
We thought about automatically changing the value of the global variable
fMRI_T0 on the basis of the user's choice of reference slice during a
previous slice-timing correction stage. However, since this change would
only apply if the same matlab session were in operation, we have decided
not change the defaults. Thus it is up to the user to remember to change
the value of fMRI_T0 (if necessary) to match the reference slice (if
slice-timing correction has been used) before creating a design matrix.
4. TA and TR
============
We originally imagined near-continuous scanning, where the TR (time
between
first slice of one scan and first slice of the next) was very close to
the
TA (time between the first and last slice of one scan). With N slices,
the time-per-slice (which the interpolation obviously needs to know) is
simply:
TR/N
To accommodate significant gaps between scans, Michael Erb allowed the
user to enter a TA smaller than the TR. Matthew Brett then pointed out a
small error in the calculation of time-per-slice, which should be (and is
now):
TA/(N-1)
Note that if your TA is a lot less than your TR, you might also want to
adjust the value of fMRI_T0 appropriately, because fMRI_T (and fMRI_T0) refer
to the number of time-bins per TR (not per TA). Thus, if you had TA=3s and
TR=4s, and wanted to reference to the middle slice in time, you would chose:
fMRI_T0 = fMRI_T * 3/4 * 1/2
Eg for the default fMRI_T=16, fMRI_T0 would be 6 (not 8).
Several people have pointed out that you could change the value of
fMRI_T to match the number of slices, N. Then the reference time-bin
could match
the reference slice (in time) exactly. This is perfectly acceptable; the
only reason we chose an fMRI_T value of 16 was that it seemed a reasonable
degree of temporal resolution with typical TRs of 3-4s: increasing fMRI_T above
this will not gain any real sensitivity (with typical TRs, given the
time constants of the HRF and temporal smoothing), but may slow down
computation.
Darren's new code (attached), which you are welcome to use, addresses
point 1 by changing to default interleaved option to a top-down odd-even
sequence, eg:
6 4 2 5 3 1
(which I repeat, may not be correct for all users - see point 1), and
addresses point 2 by asking for the reference slice in time rather than
space (which I repeat may not be relevant for interleaved sequences
- see point 2 - but does improve correspondence with the value of
fMRI_T0 - see point 3). We may include Darren's changes in the next release of
SPM.
Darren - if you are happy with this summary, please could you forward
it to the SPM mailbase? - done
Many thanks
Rik
--
---------------------------8-{)}-------------------------
DR R HENSON
Wellcome Department of Cognitive Neurology
12 Queen Square
London, WC1N 3BG
England
EMAIL: [log in to unmask]
URL: http://www.psychol.ucl.ac.uk/rik.henson/index.html
TEL1 +44 (0)20 7833 7483
TEL2 +44 (0)20 7833 7472
FAX +44 (0)20 7813 1420
---------------------------------------------------------
--
function spm_slice_timing(P, Seq, refslice, timing)
% function spm_slice_timing(P, Seq,refslice,timing)
% INPUT:
% P nimages x ? Matrix with filenames
% Seq slice acquisition order (1,2,3 = asc, desc, interl)
% 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
%
% If no input is specified the function serves as a GUI
%
% OUTPUT:
% None
%
% NMH_ACQCORRECT Correct differences in image acquisition time between slices
% 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 beginning of each TR.
%
% 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).
%
% NOTE WELL: This correction should be the first performed (i.e.,
% before orienting, motion correction, padding, smoothing, etc.).
% Additionally, it should only be performed once!
%
% Written by Darren Gitelman at Northwestern U., 1998
%
% Based (in large part) on ACQCORRECT.PRO from Geof Aquirre 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, 1999, to
% handle different sequence acquisitions, analyze format, 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 D Gitelman to make the selection of slices based on time bin.
% Note when specifying a custom order number the slices according to analyze
% convention where 1== bottom. Then specify the order as a vector reading
% left to right in which the order in the vector corresponds to the order
% of acquisition.
%_______________________________________________________________________
% @(#)spm_slice_timing.m 2.7.2 00/09/05 drg
SPMid = spm('FnBanner',mfilename,'2.7.2');
[Finter,Fgraph,CmdLine] = spm('FnUIsetup','Slice timing');
spm_help('!ContextHelp',mfilename);
nsubjects = 1;
if nargin < 1,
% Choose the images
%P = spm_get(+Inf,'*.img','Select images to acquisition correct');
% Modified by M Erb
% 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
P = [];
P = spm_get(+Inf,'*.img',...
['Select images to acquisition correct for
subject ' num2str(i)]);
eval(['P' num2str(i) ' = P;']);
end
% end of Modified by M Erb
end;
% map image files into memory
Vin = spm_vol(P);
nimgo = size(P,1);
nimg = 2^(floor(log2(nimgo))+1);
nslices = Vin(1).dim(3);
if nargin < 2,
% select fMRI acquisition sequence type
Stype = str2mat(...
'ascending (first slice=bottom)',...
'descending (first slice=top)',...
'interleaved (first slice=top)',...
'user specified');
str = 'Select sequence type';
Seq = spm_input(str,'!+1','m',Stype,[1:size(Stype,1)]);
end;
if Seq==[1],
sliceorder = [1:1:nslices];
elseif Seq==[2],
sliceorder = [nslices:-1:1];
elseif Seq==[3],
% Assumes interleaved sequences top-middle downwards
%for k=1:nslices,
% sliceorder(k) = round((nslices-k)/2 +
(rem((nslices-k),2) * (nslices - 1)/2)) + 1;
%end;
% ***** NEW CODE BY DRG *****
sliceorder = [nslices:-2:1 (nslices - 1):-2:1];
% ***** END NEW CODE *****
elseif Seq==[4],
sliceorder = [];
while length(sliceorder)~=nslices | max(sliceorder)>nslices | ...
min(sliceorder)<1 | any(diff(sort(sliceorder))~=1),
sliceorder = spm_input('Order of slices
(1=bottom)','!+0','e');
end;
end;
if nargin < 3,
% Choose reference slice (in Analyze format, slice 1 = bottom)
% Note: no checking that 1 < refslice < no.slices (default = middle slice)
% ***** MODIFIED BY DRG TO ASK FOR REFERENCE SLICE IN TIME *****
refslice = -1;
while ~any(1:nslices==refslice)
refslice = spm_input('Reference Slice in
time','!+0','e',floor(Vin(1).dim(3)/2));
end
% ***** END MODIFIED *****
end;
% ***** NEW CODE BY DRG TO PROPERLY REFERENCE REFSLICE TO SLICEORDER *****
reslice = sliceorder(refslice);
% ***** END NEW CODE *****
if nargin < 4,
%
% changed by M Erb
% factor = 1/nslices;
TR = spm_input('Interscan interval (TR) {secs}','!+1','e',3);
% TA = spm_input('Acquisition Time (TA) {secs}','!+1','e',TR);
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;
% end of changed by ME
%
else,
TR = (nslices-1)*timing(1)+timing(2);
fprintf('Your TR is %1.1f\n',TR);
factor = timing(1)/TR;
end;
spm('Pointer','Watch')
%spm('FigName','Slice timing: working',Finter,CmdLine);
% Modified by M Erb
for subj = 1:nsubjects
task=['Slice timing: working on session ' num2str(subj)];
spm('FigName',task,Finter,CmdLine);
%eval(['P = P' num2str(subj) ';']);
if nsubjects > 1, eval(['P = P' num2str(subj) ';']); end;
Vin = spm_vol(P);
nimgo = size(P,1);
nimg = 2^(floor(log2(nimgo))+1);
nslices_t= Vin(1).dim(3);
if ( nslices_t ~= nslices )
fprintf('Number of slices differ! %d %\n', nimg);
else
% end of Modified by M Erb
% create new header files
Vout = Vin;
for k=1:nimgo,
[pth,nm,xt,vr] = fileparts(deblank(Vin(k).fname));
Vout(k).fname = fullfile(pth,['a' 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)];
Vout(k) = spm_create_image(Vout(k));
end;
% 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.
for k = 1:nslices,
% Set up time acquired within slice order
shiftamount = (find(sliceorder==k) -
find(sliceorder==refslice)) * 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;
% 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
% 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])';
% 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));
% 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');
% Modified by M Erb
% endelse of "if ( nslices_t ~= nslices )"
end
% endfor of "for subj = 1:nsubjects"
end
% end of Modified by M Erb
spm('FigName','Slice timing: done',Finter,CmdLine);
spm('Pointer');
return;
--- end forwarded text
-------------------------------------------------------------------------
Darren R. Gitelman, M.D.
Cognitive Neurology and Alzheimer's Disease Center
E-mail: [log in to unmask] WWW:
http://www.brain.northwestern.edu
Voice: (312) 908-9023 Fax: (312) 908-8789
Northwestern Univ., 320 E. Superior St., Searle 11-470, Chicago, IL 60611
-------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
