topup - A tool for estimating and correcting susceptibility induced distortions

• Users guide
• Frequently asked questions

The susceptibility-induced off-resonance field

An "off-resonance" field is simply the difference between what you (or rather the scanner) think that the field is, and what it actually is. You can think about it as the field you would get if you subtracted the nominal field from the actual field. We typically think of a magnetic field in terms of Tesla, but that isn't a practical unit for these off-resonance fields that are typically on the order of micro Tesla. Instead they are typically given in Hz, where an off-resonance field of 1 Hz means that hydrogen (water) precesses 1 Hz faster than it would if there was no off-resonance field (i.e. the field was perfect).

An "off-resonance field" can be thought about as a map of the differences between what the scanner thinks the field is, and what the field actually is. For most voxels the off-resonance field will be zero, i.e. the signal from that voxel will behave according to the gradients, and as the scanner "expects" it to. But in some voxels it will be non-zero, which means that it will continously be acquiring phase in accordance with its off-resonance frequency. If for example the off-resonance is 100 Hz in a given voxel, it means that it will acquire a phase of \(2\pi\) every 0.01 seconds. The reason that there is an off-resonance field is because the head, when placed in the scanner, will slightly resist magnetisation. That means that inside the head, away from the borders or any air-filled cavities, the field will be quite flat and very slightly lower than the expected 3T (for example). On the other hand, in and around the air-filled cavities (sinuses, ear-canals etc) the field will vary quickly (spatially) and have values in excess of 100 Hz.

Distortions in EPI images

EPI images attempt to acquire all the echoes (all of k-space) it needs for an image after a single excitation. This is in contrast to many other sequences that will acquire just a single echo (line in k-space) for each excitation. This means that with EPI one can acquire a whole image in a few tens of milliseconds, and a whole volume (stack of EPI images) in a second or two. This is in contrast to for example a 3D T1-weighted image that takes several minutes to acquire. A main advantage of such a fast acquisition is the ability to follow functional changes with high temporal resolution (fMRI) and to acquire a large number of images to characterize water diffusivity (dMRI).

But it also has disadvantages, and foremost of those is the very high sensitivity to off-resonance fields. This is a consequence of the relatively long time during which the off-resonance field can translate into a phase difference. The "long time" here may seem paradoxical given that EPI is a famously fast acquisition, but the crucial time is the difference in time between the excitation (or the re-focusing pulse) and the different echoes of the k-space. For most sequences that difference is zero, i.e. the time between excitation and acquisition of the echo is exactly the same for all echoes. That means that there are no distortions (at least not in the phase-enode (PE) direction), but also that the entire sequence takes a "long" time since it involves many excitations.

With EPI, there is a single excitation followed by the acquisition of a series of echoes, each preceeded by a "PE-blip", i.e. a small phase increment. But in voxels with a non-zero off-resonance field the phase will additionally change by the product of the off-resonance field and the time between the acquisition of consecutive echoes. There is no way for the sequence/reconstruction to distinguish between this "incidental" phase accumulation and the intentional one (from the gradients), and it will lead to a displacement of signal (distortion). To make this concrete: let us assume that the PE is along the posterior<->anterior direction, and that we use positive PE-blips (meaning that the gradient causes higher field towards the top (anterior)). Hence, signal with a higher phase accumulation would indicate a location towards the top (anterior), and be located accordingly by the image reconstruction. But a voxel with a positive non-zero off-resonance field will have accumulated additional, incidental, phase and will be misplaced (distorted) towards the top (anterior). Conversely, if we used negative PE-blips (meaning that the gradient causes higher field towards the bottom (posterior)) that same voxel (with a positive off-resonance field) will now be misplaced towards the bottom.

The magnitude of the displacement (in voxels) is given by the product of the off-resonance magnitude (Hz) and the Bandwidth per voxel (Hz/voxel) in the PE-direction. The latter is a number that tells you how many Hz off-resonance that will result in signal being displaced a distance of one voxel. Typical values would be in the range 10--40 Hz/voxel. It is also the reciprocal of the total readout time (the time between the centre of the first echo and the centre of the last).

So, in conclusion, in order to know how the images are distorted one needs to know the off-resonance field, the direction and polarity of the phase-encoding and the bandwidth per voxel (or equivalently the readout time). The direction, polarity and bandwidth in the PE-direction are determined by the sequence, and are all known. The off-resonace field, by contrast, is unknown and will need to be estimated from additional measurements.

How topup estimates the off-resonance field

By acquiring two spin-echo-EPI (SE-EPI) images with opposing PE-directions one has two images with the exact same displacement of the signal at each voxel, but in opposing directions. That means that for each point \(\mathbf{x}\) in image space there is a displacement \(d\) such that the intensity that "should" be in that voxel (point) is located at a distance \(d\) and \(-d\) in the two images respectively. From that also follows follows that the intensity should be the same (save for noise) at that point \(\mathbf{x}+d\) and \(\mathbf{x}-d\) in the two images. Hence, we want to find a field \(d(\mathbf{x})\) such that for every point \(\mathbf{x}\) the displacements \(d\) and \(-d\) would lead to voxels with the same intensity in the two images. If, in addition, that field was subject to smoothness and continuity constraints it is likely that field would correspond to the "true" off-resonance field.

That is the strategy that topup uses in order to estimate the field. It is essentially a non-linear registration subject to the constraints that displacements are only along the PE-direction and of the same magnitude but opposing directions in the two input images. A pyramid approach and smoothness constraints (bending energy) are used to improve robustness and the validity of the estimates.

This means that it really estimates a displacement field, and by knowing the bandwidth it can be translated into a "proper" off-resonance field in Hz. However, that also means that it is typically not critical that the bandwidth (i.e. its reciprocal, the total readout time) is correctly specified. The displacement field is always correctly estimated, and the total readout time specification is only used to rescale it to Hz, and while it is nice to know the field in Hz it is only rarely neccessary.

One important feature of topup is that it has an internal movement model, i.e. it is able to distinguish between differences between the two images caused by subject movement and those caused by the off-resonance field. Without that all differences would be attributed to the off-resonance field, with potentially disastrous results as shown in the FAQ section. It also would not work to perform a rigid-body registration prior to estimating the off-resonance field since the two images will be very different, and the registration would not perform well. Hence, the internal movement model is a crucial part of topup.

Other ways to estimate the off-resonance field (Dual echo-time gradient echo images)

The "traditional", and possibly still the most common, method to estimate the off-resonance field is from two gradient-echo (GE) images (not EPI) acquired with different echo times (TE). In order to understand how that works we need to understand a few fundamentals of MR-physics.

Each voxel in an MR image represents a "complex number", i.e. a number that has two components referred to as the "real" and the "imaginary" part. That terminology is not necessarily helpful when it comes to gaining an intuitive understanding, and it is sufficient to think of these numbers as vectors that possess a magnitude (length) and phase (angle). In almost all MR applications the phase is irrelevant and instead an image is created where the intensity in each voxel represents the magnitude of that vector. However, for the purpose of estimating the off-resonance field using dual TE gradient-echo images it is precisely the phase that is of interest.

If the acquisition was "perfect" in every way the phase would always be zero in each voxel, and all the signal would be "real". However, that is not the case and if we first consider the short TE phase image the phase in each voxel will reflect a set of causes such as timing errors and the off-resonance field. The long TE image it will be subject to the same phase due to timing errors, but it will have been exposed to the off-resonance field for longer than the short TE image, which means that the it will have accrued phase for longer. Hence if, for a given voxel, we subtract the phase of the short TE image from the phase of the long TE image the phase from timing errors will cancel out and the difference in phase will be completely caused by the off-resonance having impacted the spins for \(\mathrm{TE}_L - \mathrm{TE}_S\) longer. Hence one can trivially calculate the off-resonace field (in Hz) as \(2 \pi (\phi_L - \phi_S) / (\mathrm{TE}_L - \mathrm{TE}_S)\).

As always in life it is not quite as easy as that. Specifically it is impossible when considering a voxel in isolation to distinguish a phase difference of \(\phi\) from one of \(\phi+2\pi n\), where \(n\) is any non-zero integer. This can only be achieved by considering a voxel in the context of neighbouring voxels and is a, susprisingly tricky, process known as "phase-unwrapping" (see prelude).

Advantages of TOPUP

Both "reverse PE-direction" and "dual echo-time" methods are able to accurately estimate the off-resonance field. However, we believe there are two definite advantages with topup. - When the field is estimated by topup it will be exactly in register with the first of the two images from which it estimates the field. Hence, if one chose the first \(b=0\) image from ones diffusion data set as the first image to pass to topup the field will be perfectly aligned with the diffusion data. It is considerably more difficult to ensure that a dual-echo time field is aligned with ones data. - It takes only ~2-3 seconds to acquire an SE-EPI volume (one in the pair with opposing PE-directions). That means that it is relatively little risk that the subject will move during the acquisition of a volume. In contrast one of the volumes in a dual-echo time pair takes tens of seconds, up to a minute, to acquire which increases the risk of subject movement. If that happens, the volume is "broken" and cannot easily be used to estimate the field.

Referencing

If you use topup in your research, please make sure that you reference at least the first of the articles listed below, and ideally the complete list. For your convenience, we provide example text (short and more details versions), which you are welcome to use in your methods description.

Brief summary text: "Data was collected with reversed phase-encode blips, resulting in pairs of images with distortions going in opposite directions. From these pairs the susceptibility-induced off-resonance field was estimated using a method similar to that described in [Andersson 2003] as implemented in FSL [Smith 2004] and the two images were combined into a single corrected one."

Andersson 2003 J.L.R. Andersson, S. Skare, J. Ashburner. How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. NeuroImage, 20(2):870-888, 2003.

Smith 2004 S.M. Smith, M. Jenkinson, M.W. Woolrich, C.F. Beckmann, T.E.J. Behrens, H. Johansen-Berg, P.R. Bannister, M. De Luca, I. Drobnjak, D.E. Flitney, R. Niazy, J. Saunders, J. Vickers, Y. Zhang, N. De Stefano, J.M. Brady, and P.M. Matthews. Advances in functional and structural MR image analysis and implementation as FSL. NeuroImage, 23(S1):208-219, 2004.

For a comparison of the performance of topup compared to "traditional" dual echo-time gradient echo fieldmaps you may also be interested in

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.