Introduction

Spatial distortions caused by susceptibility-induced field inhomogeneities are a well-known problem in echo-planar imaging (EPI) but are also starting to be investigated in 3D gradient-echo (GRE) imaging, where the push for higher resolution and lower bandwidth makes the acquisition susceptible to even small deviations from nominal B₀. This issue is particularly present in very high resolution ex vivo imaging, where the use of extremely low bandwidths can cause distortions of several (very small) pixels in the readout direction¹; and in high-resolution structural imaging², where even small artefacts can hamper the detection of longitudinal changes.

Multi-echo acquisitions are used more and more frequently in this context, in order to minimize distortions and dead time³. These datasets can be used to recover quantitative MR parameters, or combined more straightforwardly in a sum-of-squares to increase SNR. Furthermore, as odd and even echoes are acquired with an opposite readout polarity, multi-echo sequences are naturally sensitive to distortions. While these distortions between echoes currently cause misalignment and blurring, they also offer a potential avenue to estimate high-resolution fieldmaps without the need for additional data; a potential that has been untapped so far.

In contrast with TopUp⁴, which requires blip-up and blip-down images that share the same contrast, we propose Multi-Echo TopUp: a generative model of multi-echo data that jointly fits an exponential decay and the distortion map that maximizes its goodness of fit. This allows distortions to be corrected in multi-echo images without the need for additional field mapping data. This method also differs from ⁵, where the exponential model was used to synthesize reversed versions of each echo that were then matched in a TopUp manner. The proposed algorithm can be related to methods that perform model fitting and motion (or other sorts of distortions) estimation jointly^6–8.

Methods

We write the joint likelihood of the data and unknown parameters as
$$p(\mathbf{X}, \boldsymbol\Theta, \mathbf{d}) = \prod_{e=1}^{N_e} p\left(\mathbf{x}_e \mid s(\boldsymbol\Theta, t_e), \mathbf{d}\right) \times p(\boldsymbol\Theta) \times p(\mathbf{d}) ~,$$
where $$$\mathbf{X} \in \mathbb{R}^{M \times N_e}$$$ are the $$$N_e$$$ acquired echoes (with $$$M$$$ voxels), $$$\mathbf{d} \in \mathbb{R}^M$$$ is the displacement field, that gets multiplied by the parity of each echo, $$$\boldsymbol\Theta \in \mathbb{R}^{M \times 2}$$$ contains the intercept and decay parameters and $$$s$$$ is the decaying exponential function. We choose a distribution that penalizes the bending energy for $$$p(\mathbf{d})$$$, and an edge-preserving joint total-variation (JTV) distribution for $$$p(\boldsymbol\Theta)$$$. While it is common to use a Gaussian approximation for the likelihood term $$$p(\mathbf{x}_e \mid s(\boldsymbol\Theta, t_e), \mathbf{d})$$$, we found it to be too sensitive to background noise and low-SNR regions, and use a non-central Chi distribution instead. The non-central Chi likelihood is implemented in a reweighted fashion⁹, such that we are left with a nonlinear least-squares problem with respect to $$$\boldsymbol\Theta$$$ and $$$\mathbf{d}$$$. Note that in the forward model, the synthesized echo $$$s(\boldsymbol\Theta, t_e)$$$ is deformed by $$$\pm\mathbf{d}$$$ and modulated by its Jacobian determinant $$$1 \pm \nabla \mathbf{d}$$$. Our implementation alternates between optimizing the joint log-likelihood with respect to $$$\boldsymbol\Theta$$$ and $$$\mathbf{d}$$$, and uses a robust version of Gauss-Newton¹⁰ where each quadratic step is solved efficiently with a full-multi grid solver¹¹. The nonsmooth JTV prior is implemented using an iteratively reweighted scheme¹⁰.

The proposed method was tested on two ex vivo cases (whole hemisphere and cerebellum) at 7T, and two in vivo cases at 3T and 7T. Sequence parameters are provided in Table 1. All cases had fieldmaps acquired with a dual-echo GRE; and the cerebellum and in vivo cases further had reversed-polarity images acquired in an interleaved fashion² (that is, the 4 echoes are available in PNPN and NPNP orders). Our method only used one of the polarities (e.g. PNPN), whereas both versions of the first echo (i.e. P+N) were used to compute TopUp fieldmaps for comparison. We used the same TopUp parameters as in ², which were found to be optimal for fine-grained distortion correction.

Results

Figure 1 and 2 shows the result of applying all three methods to multi-echo MPRAGE data at 3T and 7T. All methods can successfully recover distortions caused by the air-tissue interface at the sphenoid sinus.
Figure 3 shows that all methods can successfully correct large distortions (up to 5 voxels) in an image of the ex vivo cerebellum.

Conclusion

We demonstrated that accurate fieldmaps can be estimated from multi-echo images without requiring a second scan with reverse polarities, and that they can be used to correct substantial distortions in low bandwidth images, both in vivo and ex vivo. This method only requires magnitude data and can be applied retrospectively to any existing multi-echo dataset.

Acknowledgements

We are grateful for the following funding sources: U01MH117023, R01AG064027, P41EB030006.