Arvind Balachandrasekaran^{1}, Merry Mani^{2}, and Mathews Jacob^{1}

Echo Planar Imaging (EPI) is widely used in many dynamic imaging studies due to its capability to provide very good temporal resolution. However, the off-resonance artifacts due to long read out result in poor correspondence with structural scans and make data interpretation difficult. Here we introduce a novel framework, where the problem of artifact correction is transformed into a recovery of image time series from undersampled measurements. We exploit the exponential structure of the signal at every pixel along with the spatial smoothness of inhomogeneity map to recover the image series. Preliminary results demonstrate the potential of the proposed method.

Echo Planar Imaging (EPI) is widely used to reduce scan time and provides good temporal resolution in many dynamic imaging studies including perfusion and functional MRI ^{1}. However, the long read-out associated with EPI makes it susceptible to field inhomogeneity induced geometric distortion artifacts. These artifacts make data interpretation difficult due to the poor correspondence of EP images with structural scans. Existing calibration-based methods estimate a field map prior to the scan and use it to correct the artifacts in the images ^{2,3}. However such methods are prone to errors due to mismatch between the actual and estimated maps due to respiration and patient motion. In this work, we propose a fast method for the calibration-free compensation of field inhomogeneity artifacts in EPI.

We adopt a time segmentation approach ^{2} and transform the original problem of 2D artifact correction into a 3D image time series recovery from undersampled measurements. We divide the acquisition window into short time segments such that the temporal evolution of phase is constant in each segment. This results in a volume with lot of missing entries as illustrated in figure (1) a. Upon recovery, the distortion-free image corresponds to the first image of the series. Under this approach, the temporal signal at every pixel location follows a single exponential model: $$ \rho[\mathbf{r},n] = \alpha(\mathbf{r}) \beta(\mathbf{r}) ^{n}; n=1,..,M $$ where $$$M$$$ is the total number of segments, $$$\alpha(\mathbf{r}) \in \mathbb{C}$$$ is the undistorted image and $$$\beta(\mathbf{r})= e^{-(\mathbf{R}_{2}^{*}(\mathbf{r})+j \omega(\mathbf{r}))T}$$$ is a spatially varying exponential function, characterized by the field map, $$$\omega(\mathbf{r})$$$, and relaxation constant, $$$\mathbf{R}_{2}^{*}(\mathbf{r})$$$. As there are two complex unknowns $$$(\alpha, \gamma = \mathbf{R}_{2}^{*} + j\omega)$$$ at every pixel location, two complex measurements are needed to estimate them. For this purpose, we acquire two EPI datasets such that the echo time (TE) of the second dataset is delayed by a few ms. The data from both the acquisitions are combined to form a new volume. This is illustrated in figure (1).

We observe that the 1D convolution of $$$\rho[\mathbf{r},n]$$$ with a two tap FIR filter yields zero. We also assume that $$$\omega(\mathbf{r})$$$ and $$$\mathbf{R}_{2}^{*}(\mathbf{r})$$$ vary smoothly across spatial locations. We thus arrive at the following annihilation relation: $$\widehat\rho[\mathbf{k},n] \otimes \widehat d [\mathbf{k},n] = 0. \tag{1}$$ where $$$\otimes$$$ denotes 3D convolution, $$$\widehat\rho[\mathbf{k},n]$$$ and $$$\widehat d [\mathbf{k},n]$$$ are the Fourier coefficients corresponding to $$$\rho[\mathbf{r},n]$$$ and the two tap filter respectively. Here the spatial bandwidth of $$$\widehat d [\mathbf{k},n]$$$ controls the smoothness of the maps, while the temporal bandwidth is two. The relation in [1] can be compactly written in terms of a multi-fold Toeplitz matrix $$$\mathcal{T}(\boldsymbol{\widehat{\rho}})$$$, which is shown to be low rank in ^{4,5}^{}. We exploit this property and propose a two step strategy to recover the image time series.

First, we estimate the null space of $$$\mathcal{T}(\boldsymbol{\widehat{\rho}})$$$ from the observed samples, as shown in figure (2). Next, we use the null space matrix $$$\mathbf{D}$$$ to recover the entire k-space volume by solving the following optimization problem: $$\min_{{\boldsymbol{\widehat{\rho}}}} \|\mathbf{S}\boldsymbol{{\widehat{\rho}}} - \mathbf{b}\|^2_2 ~ \mbox{such that}~ \mathcal{T}(\boldsymbol{\widehat{\rho}})\mathbf D = 0$$ where $$$\mathbf{S}$$$ is the sampling matrix.

We validate the proposed approach on a spherical MR phantom and two human datasets, acquired on a GE 3T scanner with a 32 channel head coil using a GRE-EPI sequence with parameters: FOV = 256mm, matrix size = 64$$$\times$$$64, slice thickness = 3.6mm, minimum TE = 30 ms, TR = 3100 ms (40 slices). For all experiments, the read out of the second dataset was delayed by $$$4 \Delta T$$$, where $$$ \Delta T$$$ = 0.636 ms. For comparison, we acquired a four shot EPI data where the distortion levels were very low.

In Figures (3), (4) and (5), we validate the proposed method on a spherical MR phantom and two human datasets. The proposed reconstructions are compared against uncorrected single shot and multi-shot EP images. In the phantom dataset, we observe that the magnetic field inhomogeneity artifacts, which cause the straight grid lines to curve in the uncorrected EP image, are corrected to a great extent in the proposed reconstruction. In the case of human datasets, we observe that the severity of the distortions are different and the proposed algorithm is able to effectively reduce them in both the datasets.

1. M. Poustchi-Amin, S. A. Mirowitz, J. J. Brown, R. C. McKinstry, and T. Li, “Principles and applications of echo-planar imaging: A review for the general radiologist,” Radiographics, vol. 21, no. 3, pp. 767–779,2001.

2. D. C. Noll, C. H. Meyer, J. M. Pauly, D. G. Nishimura, and A. Macovski, “A homogeneity correction method for magnetic resonance imaging with time-varying gradients,” IEEE Transactions on Medical Imaging, vol. 10,no. 4, pp. 629–637, 1991.

3. Bradley P Sutton, Douglas C Noll, and Jeffrey A Fessler, “Fast, iterative image reconstruction for MRI in the presence of field inhomogeneities,” IEEE Transactions on Medical Imaging, vol. 22, no. 2, pp. 178–188, 2003.

4. Arvind Balachandrasekaran, Vincent Magnotta, andMathews Jacob, “Recovery of damped exponentials using structured low rank matrix completion,” IEEE Transactions on Medical Imaging, vol. 36, no. 10, pp. 2087–2098, 2017.

5. Greg Ongie and Mathews Jacob, “Off-the-grid recovery of piecewise constant images from few Fourier samples, ”SIAM Journal on Imaging Sciences, vol. 9, no. 3, pp.1004–1041, 2016.

Figure 1: Illustration of the time segmented approach: (a) The data from the two EPI acquisitions and their corresponding uncorrected images are shown. The acquisition window of the two datasets is divided into a number of time segments. (b) By combining the time-segmented volumes of both the datasets, a k-space volume with many missing entries is formed. Data from the two EPI acquisitions lie on the yellow and green oblique planes respectively.

Figure 2: Illustration of the construction of the matrices $$$\mathcal{T}(\boldsymbol{\widehat{\rho}})$$$ and $$$\mathcal{T}_{s}(\boldsymbol{\widehat{\rho}})$$$ from the combined k-space volume $$$\boldsymbol{\widehat \rho}$$$: The 3D convolution between a filter with support $$$\Lambda$$$ and $$$\boldsymbol{\widehat \rho}$$$ results in a multi-fold Toeplitz matrix $$$\mathcal{T}(\boldsymbol{\widehat{\rho}})$$$; the valid convolutions are defined inside the red cuboid. A smaller matrix $$$\mathcal{T}_{s}(\boldsymbol{\widehat{\rho}})$$$ is constructed from $$$\mathcal{T}(\boldsymbol{\widehat{\rho}})$$$ by selecting only fully sampled rows. The null space of this smaller matrix is computed as $$$\mathbf{D}$$$.

Figure 3: Validation on phantom: The proposed reconstructions for a few slices are compared with the corresponding uncorrected single shot and multi-shot EP images in (A). We observe that the field inhomogeneity artifacts are greatly reduced in the proposed reconstruction; especially the curved grid lines in the uncorrected image appear straight in the proposed reconstruction. In (B), we also show the estimated field map and $$$\mathbf{R}_{2}^{*}$$$ maps for a particular slice.

Figure 4: Validation on human data 1: The proposed reconstructions for a few slices are compared with the corresponding uncorrected single shot and multi-shot EP images in (A). We observe that the distortions are reduced in the proposed reconstruction and are comparable to that of the multi-shot reconstruction. In (B), we also show the estimated field map and $$$\mathbf{R}_{2}^{*}$$$ maps for a particular slice.

Figure 5: Validation on human data 2: The proposed reconstructions for a few slices are compared with the corresponding uncorrected single shot and multi-shot EP images in (A). We observe that the distortions, which are more severe in this dataset, are greatly reduced in the proposed reconstruction and are comparable to that of the multi-shot reconstruction. Note that there are some differences between the proposed and the multi-shot reconstruction. This is due to the shorter echo time of the multi-shot data. In (B), we also show the estimated field map and $$$\mathbf{R}_{2}^{*}$$$ maps for a particular slice.