eddy -- a tool for correcting eddy currents and movements in diffusion data
This is a new tool to correct for eddy current-induced distortions and subject movements. It simultaneously models the effects of diffusion eddy currents and movements on the image, allowing it to work with higher b-value data than has been possible with for example eddy_correct (FSL's earlier tool for eddy current correction). As of FSL 5.0.10 (and a 5.0.9 patch) it also, optionally, performs outlier detection to identify slices where signal has been lost as a consequence of subject movement during the diffusion encoding. These slices are replaced by non-parametric predictions by the Gaussian Process that is at the heart of eddy.
Selected slice in 300 volumes of a data set from the HCP project
On the left before correction and on the right after correction for susceptibility, eddy
The diffusion signal is modelled using a Gaussian Process, which means that it makes very few assumptions about the diffusion signal (unlike parametric models such as for example the diffusion tensor). The only two assumptions are
- the signal from two acquisitions acquired with diffusion weighting along two vectors with a small angle between them is more similar than for two acquisitions with a large angle between them
the signal from two acquisitions along vectors v and -v is identical.
From these two assumptions it also follows that:
if v1 and v2 are two vectors with a "small" angle between them so that it can be assumed that the signal from the corresponding acquisitions is "similar" then v1 and -v2 are equally similar.
Because of the way the diffusion signal is modelled, and because eddy needs to be able to distinguish between signal variation caused by diffusion and that caused by eddy currents/movements it is beneficial if the data is acquired with either
A set of diffusion encoding directions that span the entire sphere and not just a half-sphere
- A blip-up-blip-down (phase encode reversed) acquisition
or with both. Note that "sampling on the whole sphere" does not imply twice as many directions. From a diffusion perspective sampling along v and -v is exactly equivalent whereas from the perspective of eddy current distortions they are different. One can therefore have two sampling schemes that both sample the diffusion evenly and equally well (and with the same total acquisition time), but where one is on the half sphere and the other on the whole sphere. To make this concrete look at the two sampling schemes below. In these plots the end of each vector is marked with an x-marker. They sample the diffusion in exactly the same way but the one on the left facilitates correcting for eddy currents while the one on the right does not. Note also that either of these schemes can easily be created from the other. To for example transform the half sphere scheme to the whole sphere scheme one just need to replace half the vectors by their negations.
Diffusion sampled on the whole sphere
Diffusion sampled on the half sphere
In order to check your own diffusion directions you can use the following Matlab commands
bvecs = load('bvecs'); % Assuming your filename is bvecs figure('position',[100 100 500 500]); plot3(bvecs(1,:),bvecs(2,:),bvecs(3,:),'*r'); axis([-1 1 -1 1 -1 1]); axis vis3d; rotate3d
The final command (rotate3d) will allow you to use the pointer to rotate the plot which is essential as from some angles a half sphere looks just like a whole sphere.
If data has neither been acquired on the whole sphere or with reversed PE-directions there is a good chance that it will still work well. In that case it can be beneficial to use the --slm=linear parameter as described in the manual.
For eddy to work well there also needs to be a minimum number of diffusion directions. The reason for this is that the concept of "close" and "distant" vectors becomes a little pointless when there are only a handful of vectors. The "minimum number" will depend on the b-value (with a larger number of directions needed for higher b-values), but it appears that the minimum is ~10-15 directions for a b-value of 1500 and ~30-40 directions for a b-value of 5000.
If it sounds like your data might be a good candidate for eddy I suggest you go on to read the manual.
If you haven't already acquired your data
If you have taken the very sensible, and unusual, step to read this documentation before acquiring your data, here are some advice.
First of all let us assume you have some maximum scan time that you cannot exceed and that this allows you to acquire N volumes. Secondly, to keep things simple, let us assume that you plan to acquire data in a single (non-zero) shell. Thirdly, you intend to run tractography on your data.
If N < 80 we recommend you acquire all your diffusion weighted images with a single PE-direction, that you acquire N unique diffusion gradients/directions optimised on the whole sphere (see above). Please note that the gradients g and -g are not unique, and that you should not acquire both. We further recommend that you acquire 2-3 b=0 volumes with an opposing PE-direction for use with topup, and that you acquire these immediately prior to the full diffusion data set. An example acquisition would be something like
PE-dir : A->P A->P P->A P->A P->A P->A P->A ...
b-value: 0 0 0 0 1500 1500 1500 ...
It can also be a good idea to intersperse additional b=0 volumes in the main data set. One in sixteen for example.
If N > 120 we suggest you consider acquiring N/2 unique diffusion gradients/directions, each acquired twice with opposing PE-directions. In this case N/2 is still > 60, which is sufficient angular sampling for most applications. The opposing PE-directions offers the option to use a different type of "interpolation" that we call "least-squares reconstruction". It works by combining the data acquired with the two PE-directions to solve the inverse problem "what might the truth look like that produce these two distorted data sets?". The method is described in detail in Andersson et al., 2003 in the reference list below. It is able to recover some of the lost resolution in areas that have been compressed by susceptibility-induced distortions.
If N is between 80 and 120 it will depend on what model you will use for your tractography, how many fibres you will model per voxel etc.
What is new in 5.0.11?
The most important new feature of eddy 5.0.11 is the ability to perform slice-to-volume motion correction. The movie below demonstrates what it can do with a problematic data set with extensive movement.
Different levels of correction of a problematic data set
Before any correction
After correction for susceptibility, eddy currents and gross subject movement
After correction for susceptibility, eddy currents, gross subject movement and signal dropout
After correction for susceptibility, eddy currents, inter- and intra-volume subject movement and signal dropout
Instead I will mention some of the changes made to the functionality associated with the --dont_sep_offs_move and --don_peas flags. These both guide the behaviour of eddy with respect to aligning the b=0 volumes with the dwi volumes. The Gaussian Process prediction maker inside eddy helps when aligning different diffusion weighted images to each other, but it is of no help when aligning dwi volumes to b=0 volumes. Hence there is need for some "post-processing" within eddy to make sure that these are aligned to each other.
There are two ways in which the b=0 and dwi volumes may end up out of alignment with each other:
- There can be a misalignment along the PE-direction caused by the difficulty to disentangle a constant eddy current-field from a subject translation along the PE-direction.
There can be an actual subject movement between the first b=0 volume and the first dwi volume. In which case the misalignment can be any combination of translations and rotations.
The former of these situations is likely to affect most data, whereas the second case is most likely to occur when there is a long gap between the first b=0 and the first dwi volume (for example when acquiring all b=0 volumes at the start) and/or the subject is uncooperative.
The first of these cases is solved by eddy by "separating the offset (constant EC field) from movement". It does so by default and only if --dont_sep_offs_move is set will eddy not do that. In previous versions this was done by attempting to work out the "offset part" using a model for gradient->EC-field. As of 5.0.11 it will also do a Mutual Information-based alignment along the PE-directon. Our tests indicate that this yields better accuracy (than the previous method) and we strongly advice against using the --dont_sep_offs_move flag.
The second of these cases is solved by performing a Mutual Information-based rigid-body registration between the mean b=0 volume and one mean volume for each dwi shell. For 5.0.11 we have changed the details of how this is done and our tests indicate that we have approximately halved the uncertainty associated with this registration. The registration is always run (even if you use the --dont_peas flag (Post Eddy Alignment of Shells)) and the rigid body parameters reported in a text file, but if the --dont_peas flag is set they will not be applied to the data.
As you will have realised by now the PEAS solution will "solve" both problems, so should one always use PEAS? I don't really have a definite answer to that question, but here are the considerations you need to make. I estimate the uncertainty of the Mutual Information-based registrations to be ~0.15 mm and ~0.15 degrees for the translations and rotations respectively.
If the first dwi volume was acquired a few seconds after the first b=0 volume and your subject was reasonably cooperative, then chances are that the actual movement between them was smaller than this uncertainty. If you do PEAS in that case you will have introduced some error (stemming from the uncertainty of PEAS) unnecessarily.
On the other hand, if you have a large number of subjects and you are not sure if maybe some of those subjects moved between the first b=0 and the first dwi volumes, then using PEAS on your entire cohort ensures that you are at least never "more wrong" than the uncertainty of the Mutual Information registration (~0.15 mm and ~0.15 degrees).
If I was tasked with analysing a data set like for example the HCP, where the first dwi volume is acquired a few seconds after the first b=0 volume and were the subjects are cooperative healthy adult, I would assume that the subject had lied still and use the --dont_peas flag. But I would also inspect the file with the PEAS parameters to ensure that there were no subjects that moved excessively in that period.
What is new in 6.0.1?
eddy 6.0.1 adds the ability to model how the susceptibility induced field changes when someone moves in the scanner. If you think of the main magnetic as a slowly flowing river (with the magnetic flux being the flow) and think of the head as an object that you insert into that river, then the disruption of the flow corresponds to the susceptibility induced field. If you translate that object across the river or upstream/downstream that river the disruption will move with the object, but its nature will be the same. This corresponds to the implicit assumption in most/all susceptibility correction methods, i.e. that as a first approximation the susceptibility field moves with the subject.
But if you now imagine rotating that object in the river (around an axis that is non-parallel with the main flux) you quickly realise that will now change the disruption it causes in a more fundamental way. The consequence of that is that a susceptibility field that you acquired/measured with the subject in one position will no longer be valid after a rotational change. The change is not big, so this is very much a "second order" correction. But in subjects/cohorts where there are lots of movement (babies, children, patients with dementia etc) it becomes important.
Static vs dynamic susceptibility correction
Before any correction
After correction for eddy currents, gross subject movement and a static susceptibility map
After correction for eddy currents, gross subject movement and a dynamic susceptibility map
As can be seen in the figure above the effects are relatively small, and mostly restricted to areas affected by susceptibility distortions in the first place (mainly frontal and temporal). It also shows that when using the dynamic susceptibility correction in eddy it as almost completely corrected.
The dynamic correction requires no additional data, just the usual static susceptibility induced fieldmap (typically calculated using ). eddy will use the data itself to calculate maps that represent the rate of change of the susceptibility field with respect to subject movement. All you need to do as a user is to include the --estimate_move_by_susceptibility flag in the eddy call. The details of this, and also some other parameters that can be set to affect how the estimation is done, can be found here.
The main reference that should be cited when using eddy is
[Andersson 2016a] Jesper L. R. Andersson and Stamatios N. Sotiropoulos. An integrated approach to correction for off-resonance effects and subject movement in diffusion MR imaging. NeuroImage, 125:1063-1078, 2016.
If you use the --repol (replace outliers) option, please also reference
[Andersson 2016b] Jesper L. R. Andersson, Mark S. Graham, Eniko Zsoldos and Stamatios N. Sotiropoulos. Incorporating outlier detection and replacement into a non-parametric framework for movement and distortion correction of diffusion MR images. NeuroImage, 141:556-572, 2016.
If you use the slice-to-volume motion model (accessed by the --mporder option) please also reference
[Andersson 2017] Jesper L. R. Andersson, Mark S. Graham, Ivana Drobnjak, Hui Zhang, Nicola Filippini and Matteo Bastiani. Towards a comprehensive framework for movement and distortion correction of diffusion MR images: Within volume movement. NeuroImage, 152:450-466, 2017.
If you use the susceptibility-by-movement correction in eddy (accessed by the option --estimate_move_by_susceptibility) please also reference
[Andersson 2018] Jesper L. R. Andersson, Mark S. Graham, Ivana Drobnjak, Hui Zhang and Jon Campbell. Susceptibility-induced distortion that varies due to motion: Correction in diffusion MR without acquiring additional data. NeuroImage, 171:277-295, 2018.
You are welcome to integrate eddy in scripts or pipelines that are subsequently made publicly available. In that case, please make it clear that users of that script should reference the paper/papers above.
Other papers of interest
For those interested in understanding the inner workings of eddy the following paper describes how it makes model-free predictions of what a diffusion weighted image should look like
[Andersson 2015] Jesper L.R. Andersson and Stamatios N. Sotiropoulos. Non-parametric representation and prediction of single- and multi-shell diffusion-weighted MRI data using Gaussian processes. NeuroImage, 122:166-176, 2015.
The "least-squares reconstruction" referred to above is described in
[Andersson 2003] Jesper L. R. Andersson, Stefan Skare and John Ashburner. How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. NeuroImage, 20:870-888, 2003.
The following paper compares the performance of eddy to the previous FSL tool eddy_correct
[Graham 2015] Mark S. Graham, Ivana Drobnjak and Hui Zhang. Realistic simulation of artefacts in diffusion MRI for validating post-processing correction techniques. NeuroImage, 125:1079-1094, 2015.
The following paper makes the case for modelling and correcting for how the susceptibility induced field changes with subject movement
[Graham 2017] M.S. Graham, I. Drobnjak, H. Zhang. Quantitative assessment of the susceptibility artefact and its interaction with motion in diffusion MRI. PLoS ONE, 12(10), 2017.