function spm_realign_ui(opt)
% Within Mode Image Realignment
%___________________________________________________________________________
%
% This routine realigns a time-series of images acquired from the same subject
% using a least squares approach and a 6 parameter (rigid body) spatial
% transformation. The first image in the list specified by the user is used
% as a reference to which all subsequent scans are realigned. The reference
% scan does not have to the the first chronologically and it may be wise to
% chose a 'representative scan' in this role.
%
% For fMRI data an additional adjustment is made to the data that removes
% a tiny amount of the movement-related confounds of these effects.
% However, it may be preferable to include the functions of the estimated
% movement parameters as confounds in the statistics part.
%
%
% Uses
% Primarily to remove movement artefact in fMRI and PET time-series (or more
% generally longitudinal studies)
%
%
% Inputs
% A series of *.img conforming to SPM data format (see 'Data Format'). The
% relative displacement of the images should be small with respect to their
% resolution. This is usually easy to ensure for functional images (e.g.
% fMRI, PET SPECT).
%
%
% Outputs
% The parameter estimation part writes out ".mat" files for each of the
% input images. The part of the routine that writes the resliced images
% uses information in these ".mat" files and writes the realigned *.img
% files to the same subdirectory prefixed with an 'r' (i.e. r*.img). The
% details of the transformation are displayed in the results window as
% plots of translation and rotation.
% A set of realignment parameters are saved for each session, named:
% realignment_params_*.txt.
%
%__________________________________________________________________________
% Refs:
%
% Friston KJ, Ashburner J, Frith CD, Poline J-B, Heather JD & Frackowiak
% RSJ (1995) Spatial registration and normalization of images Hum. Brain
% Map. 2:165-189
%
% Friston KJ, Williams SR, Howard R Frackowiak RSJ and Turner R (1995)
% Movement-related effect in fMRI time-series. Mag. Res. Med. 35:346-355
%
% W. F. Eddy, M. Fitzgerald and D. C. Noll (1996) Improved Image
% Registration by Using Fourier Interpolation. Mag. Res. Med. 36(6):923-931
%
% R. W. Cox and A. Jesmanowicz (1999) Real-Time 3D Image Registration
% for Functional MRI. Submitted to MRM (April 1999) and avaliable from:
% http://varda.biophysics.mcw.edu/~cox/index.html.
%
%_______________________________________________________________________
%
% In addition to realigning (transforming the images according to a rigid
% body model) the time series there is the option of modelling the residual
% movement related variance in an EPI (fMRI) time series that can be
% explained by a model for suceptibility-by-movement interactions. This
% is given by the "Realign & Unwarp" option.
%
% Susceptibility artefacts in EPI time series is a consequence of the
% "field disturbances" caused by the presence of an object in the field.
% Susceptibility is a property of a material, and can be thought of as
% the "resistance" put up by that material against being magnetised.
% At the interface between materials with different susceptibility
% rather severe field disturbancies will ensue. I.e. the field is
% no longer nice and homogenous. Since varying field strenght (gradients)
% is used to encode position in MRI, these disturbances cause spatial
% misplacement of RF signal, i.e. geometric distortions. Because of
% the low bandwidth in the phase encode direction these distortions
% will be appreciable mainly in that direction (typically the Ant-Post
% direction). The distortions are noticable mainly near air-tissue
% interfaces inside the scull, i.e. near the sinuses and the auditory
% canals.
%
% In these areas in particular the observed image is a severly warped
% version of reality, much like a funny mirror at a fair ground. When
% one moves in front of such a mirror ones image will distort in
% different ways and ones head may change from very elongated to
% seriously flattened. If we were to take digital snapshots of the
% reflection at these different positions it is rather obvious that
% realignment alone will not suffice to bring them into a common space.
%
% The situation is similar with EPI images, and an image collected
% for a given subject position will not be identical to that collected
% at another. We call this effect suscebtibility-by-movement interaction.
% The "Unwarp" toolbox is predicated on the assumption that the
% suscebtibility-by-movement interaction is responsible for a sizeable
% part of residual movement related variance.
%
% Assume that we know how the deformations change when the subject changes
% position (i.e. we know the derivatives of the deformations with respect
% to subject position). That means that for a given time series and a
% given set of subject movements we should be able to predict the "shape
% changes" in the object and the ensuing variance in the time series.
% It also means that, in principle, we should be able to formulate the
% inverse problem, i.e. given the observed variance (after realignment)
% and known (estimated) movements we should be able to estimate how
% deformations change with subject movement.
%
% We have made an attempt at formulating such an inverse model, and at
% solving for the "derivative fields". A deformation field can be thought
% of as little vectors at each position in space showing how that particular
% location has been deflected. A "derivative field" is then the rate of
% change of those vectors with respect to subject movement. Given these
% "derivative fields" we should be able to remove the variance caused by
% the suscebtibility-by-movement interaction. Since the underlying model
% is so restricted we would also expect experimentally induced variance
% to be preserved.
%
% Refs
%
% Background reading:
%
% Friston KJ, Williams SR, Howard R, Frackowiak RSJ and Turner R (1995)
% Movement-related effect in fMRI time-series. Mag. Res. Med. 35:346-355
%
% Jezzard P and Balaban RS (1995) Correction for geometric distortions
% in echoplanar images from B0 field variations. Magn Reson Med 34:65-73
%
% Wu DH, Lewin JS and Duerk JL (1997) Inadequacy of motion correction
% algorithms in functional MRI: Role of susceptibility-induced artefacts.
% J Magn Reson Imag 7:365-370
%
% About this particular method:
%
% Andersson JLR, Hutton C, Ashburner J, Turner R, Friston K (2001)
% Modelling geometric deformations in EPI time series. NeuroImage
% 13:903-919. doi:10.1006/nimg.2001.0746
%
%_______________________________________________________________________
%
% The Prompts Explained
%_______________________________________________________________________
%
% 'Number of subjects'
% Enter the number of subjects you wish to realign.
%
% For fMRI, it will ask you the number of sessions for each subject.
% In the coregistration step, the sessions are first realigned to
% each other, by aligning the first scan from each session to the
% first scan of the first session. Then the images within each session
% are aligned to the first image of the session.
% The parameter estimation is performed this way because it is assumed
% (rightly or not) that there may be systematic differences
% in the images between sessions.
% The adjustment step (correcting for resampling artifacts) is also
% performed completely independantly between each of the fMRI sessions.
%
% 'Images, subject # ..'
% Select the scans you wish to realign. All operations are relative
% to the first image selected.
%
% ......... Note that not all of the following prompts may be used: .........
%
% 'Which option?'
% 'Coregister only'
% Only determine the parameters required to transform each of the
% images 2..n to the same space as image 1.
% The determined parameters for image XXXX.img are saved in the
% file XXXX.mat. This is a Matlab file, containing the matrix 'M'.
% The location of an image voxel (in mm) can be determined by
% computing M(1:3,:)*[xcoord ycoord zcoord 1].
% Note that if the coregistration is performed more than once on
% the unresliced data, the starting estimates are obtained from
% parameters stored in the '.mat' files.
% Note that for PET, the coregistration is a two step process.
% First of all, the images are all realigned to the first in
% the series. A mean of these realigned images is created, and
% a second pass realignment is performed to realign all the
% images to the mean. Finally, the parameters are corrected
% for any differences estimated by registering the first image in
% the series to the mean image.
%
% 'Reslice Only'
% Reslice the specified images according to the contents of the
% previously determined parameters. The images are resliced to be
% in the same space as the first one selected. For fMRI, this is
% the first image of the first session.
%
% 'Coregister & Reslice'
% Combine the above two steps together.
%
%
% Options for reslicing:
%
% 'Create what?'
% 'All Images (1..n)'
% This reslices all the images - including the first image selected
% - which will remain in it's original position.
%
% 'Images 2..n'
% Reslices images 2..n only. Useful for if you wish to reslice
% (for example) a PET image to fit a structural MRI, without
% creating a second identical MRI volume.
%
% 'All Images + Mean Image'
% In addition to reslicing the images, it also creates a mean of the
% resliced image.
%
% 'Mean Image Only'
% Creates the mean image only.
%
%_______________________________________________________________________
%
% Defaults Options
%_______________________________________________________________________
%
% 'Registration Quality?'
% Quality versus speed trade-off. Highest quality (1) gives most
% precise results, whereas lower qualities gives faster realignment.
% The idea is that some voxels contribute little to the estimation of
% the realignment parameters. This parameter is involved in selecting
% the number of voxels that are used.
% [defaults.realign.estimate.quality]
%
% 'Allow weighting of reference image?'
% Give the option of providing a weighting image to weight each voxel
% of the reference image differently when estimating the realignment
% parameters. The weights are proportional to the inverses of the
% standard deviations.
% For example, when there is a lot of extra-brain motion - e.g., during
% speech, or when there are serious artifacts in a particular region of
% the images.
% [defaults.realign.estimate.weight]
%
%
% 'Reslice interpolation Method?'
% The method by which the images are sampled when being written in a
% different space.
% 'Nearest Neighbour'
% - Fastest, but not normally recommended.
% 'Bilinear Interpolation'
% - OK for PET, or realigned fMRI.
%
% 'B-spline Interpolation'
% - Better quality (but slower) interpolation, especially
% with higher degree splines. Don't use B-splines when
% there is any region of NaN or Inf in the images.
% 'Fourier space Interpolation'
% - Rigid body rotations are executed as a series of shears,
% which are performed in Fourier space (Eddy et. al. 1996).
% Unfortunately, this method can only be applied to images
% with cubic voxels (since zooms can not be done by
% convolution in Fourier space).
% [defaults.realign.write.interp]
%
% These are typically:
% 'No wrapping' - for PET or images that have already
% been spatially transformed.
% 'Wrap in Y' - for (un-resliced) MRI where phase encoding
% is in the Y direction (voxel space).
% [defaults.realign.write.wrap]
%
% 'Mask images?'
% Because of subject motion, different images are likely to have different
% patterns of zeros from where it was not possible to sample data.
% With masking enabled, the program searches through the whole time series
% looking for voxels which need to be sampled from outside the original
% images. Where this occurs, that voxel is set to zero for the whole set
% of images (unless the image format can represent NaN, in which case
% NaNs are used where possible). This is in order to avoid artifactual
% movement-related variance the realigned images.
% [defaults.realign.write.mask]
%
% Other Defaults are useful for making spm_realign_ui.m more flexible:
% [defaults.realign.estimate.interp] - interpolation method for parameter
% estimation.
% [defaults.realign.estimate.wrap] - wrapping used for parameter estimation.
% This should probably also influence the way that smoothing is done, but
% the smoothing code has not yet been modified to include wrapping.
%
%__________________________________________________________________________
%
% The `.mat' files.
%
% This simply contains a 4x4 affine transformation matrix in a variable `M'.
% These files are normally generated by the `realignment' and
% `coregistration' modules. What these matrixes contain is a mapping from
% the voxel coordinates (x0,y0,z0) (where the first voxel is at coordinate
% (1,1,1)), to coordinates in millimeters (x1,y1,z1). By default, the
% the new coordinate system is derived from the `origin' and `vox' fields
% of the image header.
%
% x1 = M(1,1)*x0 + M(1,2)*y0 + M(1,3)*z0 + M(1,4)
% y1 = M(2,1)*x0 + M(2,2)*y0 + M(2,3)*z0 + M(2,4)
% z1 = M(3,1)*x0 + M(3,2)*y0 + M(3,3)*z0 + M(3,4)
%
% Assuming that image1 has a transformation matrix M1, and image2 has a
% transformation matrix M2, the mapping from image1 to image2 is: M2\M1
% (ie. from the coordinate system of image1 into millimeters, followed
% by a mapping from millimeters into the space of image2).
%
% These `.mat' files allow several realignment or coregistration steps to be
% combined into a single operation (without the necessity of resampling the
% images several times). The `.mat' files are also used by the spatial
% normalisation module.
%__________________________________________________________________________
%
% The *uw.mat file.
%
% The *uw.mat file contains information about input to the Unwarp module
% (i.e. how it was run) and about the result. It contains all the information
% needed to unwarp the appurtenant time-series, and to assess exactly how the
% estimation was done.
%__________________________________________________________________________
% @(#)spm_realign_ui.m 2.13 John Ashburner - with input from Oliver Josephs 03/01/24
% and Jesper Andersson
global defaults
if nargin==0 | strcmp(lower(opt),'ui'),
if get(gcbo,'Value') < 2
run_ui(defaults.realign, defaults.modality);
else
run_ui(defaults.realign, defaults.modality, defaults.unwarp);
end
elseif nargin>0 & strcmp(lower(opt),'defaults'),
defaults.realign = get_defs(defaults.realign);
elseif nargin>0 & strcmp(lower(opt),'unwarpdefaults'),
defaults.unwarp = get_unwarp_defs(defaults.unwarp);
end;
return;
function run_ui(defs, modality, unwarp)
% User interface.
%_______________________________________________________________________
SPMid = spm('FnBanner',mfilename,'2.10');
[Finter,Fgraph,CmdLine] = spm('FnUIsetup','Realign');
spm_help('!ContextHelp',mfilename);
n = spm_input('Num subjects', '+1', 'e', 1);
if n<1, spm_figure('Clear','Interactive'); return; end
%
% Possibility to include measured field-map temporarily removed.
%
pcpm = 0;
%if strcmp(lower(modality),'fmri') & nargin > 2
% pcpm = spm_input('Include pre-calculated phase maps?','+1','m',...
% [' Yes| No'],[1 0],0);
%end
P = cell(n,1);
for i = 1:n,
if strcmp(lower(modality),'fmri'),
ns = spm_input(['Num sessions, subj ' num2str(i)], '+1', 'e', 1);
pp = cell(1,ns);
for s=1:ns,
p = '';
pm = '';
while size(p,1)<1,
p = spm_get(Inf,'IMAGE',...
['Images, subj ' num2str(i) ', sess' num2str(s)]);
if pcpm == 1
pm = spm_get(1,'IMAGE',...
['Phasemap, subj ' num2str(i) ', sess' num2str(s)]);
end
end;
pp{s} = p;
ppm{s} = pm;
end;
P{i} = pp;
pP{i} = ppm;
else,
p = cell(1,1);
p{1} = '';
while size(p{1},1)<1,
p{1} = spm_get(Inf,'IMAGE',...
['Images, subj ' num2str(i)]);
end;
P{i} = p;
end;
end;
if strcmp(lower(modality),'pet'),
FlagsC = struct('quality',defs.estimate.quality,'fwhm',7,'rtm',1);
else,
FlagsC = struct('quality',defs.estimate.quality,'fwhm',5,'rtm',0);
end;
if nargin < 3
WhchPtn = spm_input('Which option?', '+1', 'm',...
'Coregister only|Reslice Only|Coregister & Reslice',...
[1 2 3],3);
else
WhchPtn = 3;
end
PW = '';
if (WhchPtn == 1 | WhchPtn == 3) & defs.estimate.weight,
if spm_input(...
['Weight the reference image(s)?'],...
'+1', 'm',...
['Dont weight registration|'...
'Weight registration'], [0 1], 1);
PW = spm_get(n,'IMAGE', 'Weight images for each subj');
end;
end;
if nargin > 2
foe = spm_input('Model field changes w.r.t.','+1','m',...
['Pitch & Roll|All|Customise'],...
[1:2 0],1);
if foe == 1
foe = [4 5];
elseif foe == 2
foe = 1:6;
else
foe = spm_input('First order effects','+0','e',...
'4 5',Inf)';
end
if unwarp.estimate.soe == 2
tmp = spm_input('Second order effects','+1','m',...
['None|All first order effects|Customise'],...
[1:3],1);
if tmp == 1
soe = [];
else
cnt = 1;
for i=1:size(foe,2)
for j=i:size(foe,2)
soe(cnt,1) = foe(i);
soe(cnt,2) = foe(j);
cnt = cnt+1;
end
end
if tmp == 3
string = '[';
for i=1:cnt-1
string = [string sprintf('%d %d; ',soe(i,1),soe(i,2))];
end
string = [string(1:end-2) ']'];
soe = spm_input('Second order effects','+0','e',...
string,[Inf 2]);
end
end
else
soe = [];
end
end
% Reslicing options
%-----------------------------------------------------------------------
if WhchPtn == 2 | WhchPtn == 3,
FlagsR = struct('interp',defs.write.interp,...
'wrap',defs.write.wrap,...
'mask',defs.write.mask,...
'which',2,'mean',1);
if strcmp(lower(modality),'pet'), FlagsR.wrap = [0 0 0]; end;
if nargin < 3
p = spm_input('Create what?','+1','m',...
[' All Images (1..n)| Images 2..n|'...
' All Images + Mean Image| Mean Image Only'],...
[1 2 3 4],3);
if p==1, FlagsR.which = 2; FlagsR.mean = 0; end
if p==2, FlagsR.which = 1; FlagsR.mean = 0; end
if p==3, FlagsR.which = 2; FlagsR.mean = 1; end
if p==4, FlagsR.which = 0; FlagsR.mean = 1; end
else
p = spm_input('Create what?','+1','m',...
[' All Images (1..n)| All Images + Mean Image'],[1 3],3);
if p==1, FlagsR.which = 2; FlagsR.mean = 0; end
if p==3, FlagsR.which = 2; FlagsR.mean = 1; end
end
end;
%
% If doing realignment only
%
if nargin < 3
spm('Pointer','Watch');
for i = 1:n
if WhchPtn==1 | WhchPtn==3,
spm('FigName',['Realigning subj ' num2str(i)],Finter,CmdLine);
flagsC = FlagsC;
if ~isempty(PW), flagsC.PW = deblank(PW(i,:)); end;
spm_realign(P{i},flagsC);
end
if WhchPtn==2 | WhchPtn==3,
spm('FigName',['Reslicing subj ' num2str(i)],Finter,CmdLine);
spm_reslice(P{i},FlagsR);
end;
end;
spm('FigName','Realign: done',Finter,CmdLine);
spm('Pointer');
%
% If doing realignment and unwarping
%
else
uwe_flags = struct('order', unwarp.estimate.basfcn,...
'sfield', [],...
'regorder', unwarp.estimate.regorder,...
'lambda', unwarp.estimate.regwgt,...
'jm', unwarp.estimate.jm,...
'fot', foe,...
'sot', soe,...
'fwhm', unwarp.estimate.fwhm,...
'rem', unwarp.estimate.rem,...
'noi', unwarp.estimate.noi,...
'exp_round', unwarp.estimate.expround);
uwr_flags = struct('interp', defs.write.interp,...
'wrap', defs.write.wrap,...
'mask', defs.write.mask,...
'which', FlagsR.which,...
'mean', FlagsR.mean);
if unwarp.estimate.jm == 1
uwr_flags.udc = 2;
else
uwr_flags.udc = 1;
end
for i = 1:n
spm('Pointer','Watch');
spm('FigName',['Realigning subj ' num2str(i)],Finter,CmdLine);
flagsC = FlagsC;
if ~isempty(PW), flagsC.PW = deblank(PW(i,:)); end;
spm_realign(P{i},flagsC);
spm('Pointer');
clear ads
tmpP = spm_vol(P{i}{1}(1,:));
uwe_flags.M = tmpP.mat;
for j=1:length(P{i})
if pcpm > 0
uwe_flags.sfield = pP{i}{j};
end
ds = spm_uw_estimate(P{i}{j},uwe_flags);
ads(j) = ds;
[path,name,ext,ver] = fileparts(P{i}{j}(1,:));
pefile = fullfile(path,[name '_uw.mat']);
save(pefile,'ds');
end
spm_uw_apply(ads,uwr_flags);
end
end
return;
%_______________________________________________________________________
%_______________________________________________________________________
function defs = get_defs(defs)
tmp2 = [1.00 0.90 0.75 0.50 0.25 0.10 0.05 0.01 0.005 0.001];
tmp = find(defs.estimate.quality == tmp2);
if isempty(tmp) tmp = 1; end
defs.estimate.quality = spm_input('Registration Quality?','+1','m',...
['Quality 1.00 (slowest/most accurate) |Quality 0.90|' ...
'Quality 0.75|Quality 0.50|Quality 0.25|Quality 0.10|' ...
'Quality 0.05|Quality 0.01|' ...
'Quality 0.005|Quality 0.001 (fastest/poorest)'],tmp2, tmp);
tmp = 1;
if defs.estimate.weight == 1, tmp = 2; end;
defs.estimate.weight = spm_input(...
['Allow weighting of reference image?'],...
'+1', 'm',...
['Dont allow weighting|'...
'Allow weighting'], [0 1], tmp);
tmp2 = [0 1 2 3 4 5 6 7 Inf];
tmp = find(defs.write.interp == tmp2);
if ~finite(defs.write.interp), tmp = 9; end;
if isempty(tmp), tmp = 2; end;
defs.write.interp = spm_input('Reslice interpolation method?','+1','m',...
['Nearest Neighbour|Trilinear|2nd Degree B-Spline|'...
'3rd Degree B-Spline|4th Degree B-Spline|5th Degree B-Spline|'...
'6th Degree B-Spline|7th Degree B-Spline|Fourier Interpolation'],...
tmp2,tmp);
wraps = [0 0 0 ; 1 0 0; 0 1 0; 1 1 0; 0 0 1; 1 0 1; 0 1 1; 1 1 1];
t = find(all(repmat(defs.write.wrap(:)',8,1) == wraps, 2));
if isempty(t), t = 1; end;
p = spm_input('Way to wrap images?','+1','m',...
['No wrap|Wrap X|Wrap Y|Wrap X & Y|Wrap Z|Wrap X & Z|Wrap Y & Z|Wrap X, Y & Z'],...
[1 2 3 4 5 6 7 8], t);
defs.write.wrap = wraps(p,:);
defs.estimate.wrap = defs.write.wrap;
tmp = 1;
if ~defs.write.mask, tmp = 2; end;
defs.write.mask = spm_input(['Mask images?'], '+1', 'm',...
' Mask images|Dont mask images', [1 0], tmp);
return;
%_______________________________________________________________________
%_______________________________________________________________________
function defs = get_unwarp_defs(defs)
orders = [6 6; 8 8; 10 10; 12 12];
lambdas = [1e4 1e5 1e6];
defs.estimate.fwhm = spm_input('Filter width (mm)','+1','e',...
'4',1);
defs.estimate.basfcn = spm_input('No. of basis functions','+1','m',...
['6x6x*|8x8x*|10x10x*|12x12x*'],...
[1:4],3);
defs.estimate.basfcn = orders(defs.estimate.basfcn,:);
defs.estimate.regwgt = spm_input('Amount of regularisation','+1','m',...
['A little|Medium|A lot'],...
[1:3],2);
defs.estimate.regwgt = lambdas(defs.estimate.regwgt);
defs.estimate.rem = spm_input('Re-estimation of movement parameters?','+1','m',...
['Yes|No'],[1 0],1);
defs.estimate.jm = spm_input('Include Jacobian intensity modulation?','+1','m',...
['Yes|No'],[1 0],2);
defs.estimate.soe = spm_input('Second order effects','+1','m',...
['None|Customise'],...
[1:2],1);
return