Fan Lam^{1,2} and Bradley P. Sutton^{1,2}

Intravoxel macroscopic B_{0} field inhomogeneity is common source for image artifacts, and can cause significant errors in quantitative analysis of the imaging data, particularly for gradient-echo-based acquisitions. This work presents a new method for intravoxel B_{0} inhomogeneity corrected reconstruction. Specifically, we derived a signal model that incorporates intravoxel B_{0} variations, and introduced a low-rank approximation to the encoding operator. This approximation allows for very efficient computation of the forward model, which can be integrated into any regularized reconstruction formulation. Effective correction for multi-gradient-echo data were demonstrated. We expect the proposed method to be useful for a range of applications involving T_{2}^{*} contrast.

A general signal equation for multi-gradient-echo acquisition can be written as follows (ignoring the phase effects during each short readout period)

$$d_c(\mathbf{k},t)=\int_{-\infty}^{\infty}s_c(\mathbf{r})\rho(\mathbf{r},t)e^{i2\pi\Delta f(\mathbf{r})t}e^{-i2\pi\mathbf{kr}}d\mathbf{r} \qquad\qquad\quad\quad [1]$$

where $$$\rho(\mathbf{r},t)$$$ is the multi-echo images with $$$t$$$ the echo times, $$$\Delta f(\mathbf{r})$$$ the B_{0} inhomogeneity, $$$s_c(\mathbf{r})$$$ the coil sensitivity, $$$\mathbf{r}$$$ and $$$\mathbf{k}$$$ vectors denoting the spatial and k-space coordinates. Discretizing Eq. [1] using rectangular voxel basis functions, we can obtain^{4}

$$d_c(\mathbf{k}_m,t)=\sum_ns_{c,n}\boldsymbol{\rho}_n(t)e^{-i2\pi\mathbf{k}_m\mathbf{r}_n}\mathbf{W}_{mn}(t) \qquad\qquad\qquad\quad [2]$$

where $$$\{s_{c,n}\}$$$ and $$$\{\boldsymbol{\rho}_n(t)\}$$$ represent the sensitivity maps and the unknown images, respectively. The $$$M\times N$$$ matrix $$$\mathbf{W}(t)$$$ encodes the effects of B_{0} inhomogeneity (echo dependent). $$$M$$$ is the total number of k-space samples each echo and $$$N$$$ is the number of voxels. With a locally linear assumption on $$$\Delta f(\mathbf{r})$$$, $$$\mathbf{W}(t)$$$ can be simplified to $$$\mathbf{W}_{mn}(t)=e^{i2\pi f_nt}\text{sinc}[(\mathbf{k}_m-\mathbf{g}_nt)\Delta\mathbf{r}]$$$ with $$$\mathbf{g}_n$$$ being 3x1 vectors capturing the within-voxel field gradients along all three spatial directions.

It can be seen that 1) storing the matrix $$$\mathbf{W}$$$ can be memory prohibitive and 2) directly evaluating Eq. [2] is rather computationally intenstive, especially for high-resolution 3D reconstruction, because the coupling between $$$\mathbf{W}$$$ and $$$e^{-i2\pi\mathbf{k}_m\mathbf{r}_n}$$$ prevents direct application of FFT. However, we recognize that $$$\mathbf{W}$$$ can be accurately approximated by a low-rank factorization (illustrated in Fig. 1), i.e., $$$\mathbf{W}(t)\approx\sum_{l=1}^L\mathbf{u}_{l}(t)\mathbf{v}_{l}^T(t)$$$ where $$$\mathbf{u}$$$, $$$\mathbf{v}$$$ are column vectors and $$$L\ll M,N$$$. Accordingly, Eq. [2] can be rewritten in the following matrix-vector multiplication form

$$\mathbf{d}_{c,t}=\sum_{l=1}^L\mathbf{u}_l(t)\odot\{\mathbf{F}_{\Omega}(\mathbf{v}_l(t)\odot\mathbf{S}_c\boldsymbol{\rho}_{t})\}=\mathbf{E}_{c,t}\boldsymbol{\rho}_{t} \qquad\qquad\quad [3]$$

with $$$\mathbf{d}_{t}$$$ containing data for the echo image $$$\boldsymbol{\rho}_{t}$$$, $$$\mathbf{F}_{\Omega}$$$ denoting the Fourier transform operator with sampling mask $$$\Omega$$$, $$$\odot$$$ the point-wise multiplication and $$$\mathbf{S}_c$$$ the sensitivity modulation. This low-rank decomposed encoding operator in Eq. [3] can be readily incorporated into any reconstruction formulation and calculated efficiently (with L FFTs). Therefore, the final intravoxel B_{0} inhomogeneity corrected multi-echo images are jointly reconstructed by solving the following regularized least-squares problem

$$\{\hat{\boldsymbol{\rho}}_{t}\}=\arg\underset{\{\boldsymbol{\rho}_{t}\}}{\min}\sum_{t}\sum_{c}\left\Vert\mathbf{d}_{c,t}-\mathbf{E}_{c,t}\boldsymbol{\rho}_{t}\right\Vert_2^2+\lambda R(\{\boldsymbol{\rho}_{t}\}),$$

for which different fast algorithms can be used for different choices of $$$R(\cdot)$$$. Note that for each t, the factorization $$$\{\mathbf{u}_l\}$$$ and $$$\{\mathbf{v}_l\}$$$ can be computed using only sparse samples in $$$\mathbf{W}$$$ (due to low-rankness), thus avoiding explicitly storing the huge $$$\mathbf{W}$$$ in memory, making efficient implementations possible^{7,8}. Table 1 further highlights the memory and computational complexity reduction offered by the proposed method.

In vivo experiments were performed to evaluate the proposed method. First, a high-resolution (HR) multi-echo 3D GRE data (0.9x0.9x1.5mm^{3}) was acquired on a 3T Siemens scanner using a 20-channel head coil. Lower-resolution (LR) data was generated by a k-space truncation (to a 2x2x3mm^{3} voxel size). Figure 2 compares a reconstruction from the HR data (phase removed at the original resolution followed by truncating the k-space data), a Fourier reconstruction from the LR data, and the proposed reconstruction from the LR data with $$$L=8$$$. As can be seen, the HR data has minimal intravoxel dephasing due to the small voxel size, while the LR Fourier reconstruction has apparent artifacts in regions with strong susceptibility variations. The proposed method significantly reduced these artifacts (Fig. 2c), producing similar images as the truncated magnitude HR data. The R_{2}* maps shown in Fig. 3 further demonstrate the utility of the proposed method in improving the accuracy of quantitative analysis for multi-gradient-echo data.

1. Ordidge et al., Assessment of relative brain iron concentrations using T_{2}‐weighted and T_{2}*‐weighted MRI at 3 Tesla. Magn. Reson. Med., 1994;32:335-341.

2. Fernandez‐Seara and Wehrli, Postprocessing technique to correct for background gradients in image‐based R_{2}^{*} measurements. Magn. Reson. Med., 2000;44:358-366.

3. Hernando et al., R_{2}^{* }mapping in the presence of macroscopic B_{0} field variations. Magn. Reson. Med., 2012;68:830-840.

4. Yablonskiy et al., Voxel spread function method for correction of magnetic field inhomogeneity effects in quantitative gradient-echo-based MRI. Magn. Reson. Med., 2013;70:1283-1292.

5. Sutton et al., Fast, iterative image reconstruction for MRI in the presence of field inhomogeneities. IEEE Trans. Med. Imaging, 2003;22:178-188.

6. Ngo and Sutton, R_{2}^{*} mapping for robust brain function detection in the presence of magnetic field inhomogeneity. IEEE-EMBC, 2014.

7. Liang, Spatiotemporal imaging with partially separable functions. IEEE-ISBI, 2007.

8. Haldar and Hernando, Rank-constrained solutions to linear matrix equations using power factorization, IEEE Signal Process. Lett., 2009;16:584-587.

Figure 1: An illustration of the within-voxel B_{0} field variations and the low rankness of the matrix $$$\mathbf{W}$$$ in Eq. [2] that incorporates such effects. Image (a) shows the z-direction intravoxel field gradients (in Hz/mm) for the multi-echo data. Image (b) shows the singular value decay for $$$\mathbf{W}$$$ with $$$\mathbf{W}_{mn}(t)=\text{sinc}[(\mathbf{k}_m-\mathbf{g}_nt)\Delta\mathbf{r}]$$$ ($$$\mathbf{g}_n$$$ the intravoxel field gradients). The dash line indicates the threshold beyond which the rank truncation error is less than 1% ($$$L=8$$$ in this case).

Table 1: A comparison of memory usage and computational complexity between a direction evaluation of Eq. [2] with explicit formation of $$$\mathbf{W}$$$ and the proposed low-rank operator. With a small $$$L$$$, significant reduction can be achieved. Specifically, for a 128x128x32 matrix size, i.e., $$$N=M=$$$128x128x32, the memory usage can be reduced from 2 **TB** (double precision) to 64 **MB** with a rank 8 approximation, and the computational complexity can be reduced from $$$O(10^{11})$$$ to approximately $$$O(10^{8})$$$.

Figure 2: Reconstructions from the 3D GRE data for different slices across the 3D volume (7th echo). Images were produced by Fourier reconstruction of the high-resolution data (a), Fourier reconstruction of the low-resolution data (b), and the proposed reconstruction of the low-resolution data (c). The images in (a) were obtained by taking magnitude of the high-resolution reconstruction followed by a k-space truncation. Significant reduction of intravoxel dephasing induced distortion was achieved by the proposed method.

Figure 3: R_{2}^{*} maps estimated from (a) the Fourier reconstruction of the low-resolution data (corresponding to Fig. 2b) and (b) the proposed reconstruction (corresponding to Fig. 2c). Substantial reduction of intravoxel dephasing induced estimation error can be observed for the proposed reconstruction, as indicated by the red arrows.