JiscMail Logo
Email discussion lists for the UK Education and Research communities

Help for SPM Archives


SPM Archives

SPM Archives


SPM@JISCMAIL.AC.UK


View:

Message:

[

First

|

Previous

|

Next

|

Last

]

By Topic:

[

First

|

Previous

|

Next

|

Last

]

By Author:

[

First

|

Previous

|

Next

|

Last

]

Font:

Proportional Font

LISTSERV Archives

LISTSERV Archives

SPM Home

SPM Home

SPM  2000

SPM 2000

Options

Subscribe or Unsubscribe

Subscribe or Unsubscribe

Log In

Log In

Get Password

Get Password

Subject:

Fwd: SLICE-TIMING - FINAL WORD?

From:

Darren Gitelman <[log in to unmask]>

Reply-To:

Darren Gitelman <[log in to unmask]>

Date:

Wed, 13 Sep 2000 23:49:07 -0500

Content-Type:

text/plain

Parts/Attachments:

Parts/Attachments

text/plain (578 lines)

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
-------------------------------------------------------------------------


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Top of Message | Previous Page | Permalink

JiscMail Tools


RSS Feeds and Sharing


Advanced Options


Archives

September 2022
August 2022
July 2022
June 2022
May 2022
April 2022
March 2022
February 2022
January 2022
December 2021
November 2021
October 2021
September 2021
August 2021
July 2021
June 2021
May 2021
April 2021
March 2021
February 2021
January 2021
December 2020
November 2020
October 2020
September 2020
August 2020
July 2020
June 2020
May 2020
April 2020
March 2020
February 2020
January 2020
December 2019
November 2019
October 2019
September 2019
August 2019
July 2019
June 2019
May 2019
April 2019
March 2019
February 2019
January 2019
December 2018
November 2018
October 2018
September 2018
August 2018
July 2018
June 2018
May 2018
April 2018
March 2018
February 2018
January 2018
December 2017
November 2017
October 2017
September 2017
August 2017
July 2017
June 2017
May 2017
April 2017
March 2017
February 2017
January 2017
December 2016
November 2016
October 2016
September 2016
August 2016
July 2016
June 2016
May 2016
April 2016
March 2016
February 2016
January 2016
December 2015
November 2015
October 2015
September 2015
August 2015
July 2015
June 2015
May 2015
April 2015
March 2015
February 2015
January 2015
December 2014
November 2014
October 2014
September 2014
August 2014
July 2014
June 2014
May 2014
April 2014
March 2014
February 2014
January 2014
December 2013
November 2013
October 2013
September 2013
August 2013
July 2013
June 2013
May 2013
April 2013
March 2013
February 2013
January 2013
December 2012
November 2012
October 2012
September 2012
August 2012
July 2012
June 2012
May 2012
April 2012
March 2012
February 2012
January 2012
December 2011
November 2011
October 2011
September 2011
August 2011
July 2011
June 2011
May 2011
April 2011
March 2011
February 2011
January 2011
December 2010
November 2010
October 2010
September 2010
August 2010
July 2010
June 2010
May 2010
April 2010
March 2010
February 2010
January 2010
December 2009
November 2009
October 2009
September 2009
August 2009
July 2009
June 2009
May 2009
April 2009
March 2009
February 2009
January 2009
December 2008
November 2008
October 2008
September 2008
August 2008
July 2008
June 2008
May 2008
April 2008
March 2008
February 2008
January 2008
December 2007
November 2007
October 2007
September 2007
August 2007
July 2007
June 2007
May 2007
April 2007
March 2007
February 2007
January 2007
2006
2005
2004
2003
2002
2001
2000
1999
1998


JiscMail is a Jisc service.

View our service policies at https://www.jiscmail.ac.uk/policyandsecurity/ and Jisc's privacy policy at https://www.jisc.ac.uk/website/privacy-notice

For help and support help@jisc.ac.uk

Secured by F-Secure Anti-Virus CataList Email List Search Powered by the LISTSERV Email List Manager