Autocalibrating and Calibrationless Parallel Magnetic Resonance Imaging as a Bilinear Inverse Problem

TL;DR

This paper discusses recent approaches to jointly estimate unknown coil sensitivities and images in MRI, formulating it as a bilinear inverse problem to improve undersampled image reconstruction.

Contribution

It introduces and reviews methods for solving the bilinear inverse problem in MRI, enabling calibrationless and auto-calibrating image reconstruction.

Findings

Enhanced image quality from undersampled data

Reduced need for separate calibration scans

Improved robustness to subject movement

Abstract

Modern reconstruction methods for magnetic resonance imaging (MRI) exploit the spatially varying sensitivity profiles of receive-coil arrays as additional source of information. This allows to reduce the number of time-consuming Fourier-encoding steps by undersampling. The receive sensitivities are a priori unknown and influenced by geometry and electric properties of the (moving) subject. For optimal results, they need to be estimated jointly with the image from the same undersampled measurement data. Formulated as an inverse problem, this leads to a bilinear reconstruction problem related to multi-channel blind deconvolution. In this work, we will discuss some recently developed approaches for the solution of this problem.

Figures

Equations

$$y_j(\vec{k})={\int }_{\Omega }\text{d}\vec{r}\rho (\vec{r})c_j(\vec{r})e^{-i2\pi \vec{k}\vec{r}}$$

$$K_{st}(\vec{u},\vec{v})={\int }_{\Omega }\text{d}\vec{r}{c}_s(\vec{r})\overline{c_t(\vec{r})}{e}^{-i2\pi \left(\vec{u}-\vec{v}\right)\vec{r}}.$$

$$\underset{x}{\mathrm{argmin}}\sum _j\Vert P\mathcal{F}{S}_{j}x-y_j{\Vert }_2^2+R(x)$$

$$\sum _r{R}_r^H{V}_{\perp }{V}_{\perp }^H{R}_r\hat{y}$$

$$\underset{\mathcal{W}}{\underbrace{M^{-1}\sum _r{R}_r^H{V}_{\parallel }{V}_{\parallel }^H{R}_r}}y$$

$$\underset{x^i}{\mathrm{argmin}}\sum _j\Vert P\mathcal{F}\sum _i{S}_j^i{x}^i-y_j{\Vert }_2^2+\sum _{i}R(x^i)$$

$$\underset{\hat{y}}{\mathrm{argmin}}\Vert y-P\hat{y}\Vert \text{with}\mathrm{rank}\mathcal{C}\left\{\hat{y}\right\}\le L$$

$$\underset{x,c_j}{\mathrm{argmin}}$$

$$\frac{\alpha }{2}(\sum _j\Vert W{c}_j{\Vert }_2^2+\Vert x{\Vert }_2^2)$$

$$\underset{x,c_j}{\mathrm{argmin}}$$

$$\frac{\alpha }{2}(\sum _{i,j}\Vert W{c}_j^i{\Vert }_2^2+\sum _i\Vert x^i{\Vert }_2^2)$$

Full text

I Parallel Imaging

The multi-channel signal $y_{j}$ in parallel MRI is given by the Fourier transform of magnetization image $\rho$ weighted by the sensitivities of the receive coils $c_{j}$:

[TABLE]

Here, $\vec{k}$ is the sampling trajectory in the Fourier domain (k-space), and $\Omega$ is a compactly supported region inside $\mathcal{R}^{2}$ for two-dimensional or $\mathcal{R}^{3}$ for three-dimensional imaging experiments. The extend of $\Omega$ (field-of-view) in the image domain defines the bandwidth of the signal in k-space. If k-space is sampled according to the corresponding Nyquist rate it is referred to as fully sampled. As the signal equation is valid only for short durations, the acquisition of all required k-space samples has to be split into a repeated series of individual measurements, which makes the measurement process very time consuming. The main objective in parallel MRI is to reconstruct diagnostic images while reducing the number of samples as much as possible. Fig. 1 shows various common sampling schemes and Fig. 2 shows individual coil images and sensitivities for a modern 32-channel receiver-coil array.

The signals for all channels are highly correlated. In fact, it can be shown that the space of ideal signals $y_{j}(\vec{k})$ is a subspace of all band-limited k-space signals. Parallel imaging exploits this to reduce the number of samples below the Nyquist limit. The subspace of ideal signals can be characterized as a Reproducing Kernel Hilbert Space with a matrix-valued reproducing kernel given by:

[TABLE]

For known sensitivities $c_{j}$, the reproducing kernel can be computed and standard approximation theory then yields a complete sampling theory with optimal formulas for interpolation and error bounds [1].

More commonly, image reconstruction is formulated as a linear inverse problem and a regularized solution is defined as the minimizer of the following optimization problem [2, 3]:

[TABLE]

Here, $S_{j}$ is the multiplication with the sensitivity $c_{j}$, $\mathcal{F}$ the Fourier transform, and $P$ the sampling operator, and $R(x)$ a regularization term. Most state-of-the art algorithms for MRI are based on numerical optimization of the discretized functional and often also combine parallel imaging with concepts from compressed sensing to further reduce the number of samples [4, 5].

II Autocalibrating Parallel Imaging

In practice, the coil sensitivities $c_{j}$ are not known. Many techniques have been developed for autocalibration of the sensitivities from a fully-sampled (calbration) region $y_{calib}$ in the center of k-space parallel imaging. Early schemes estimated the coil sensitivities directly from low-resolution images obtained as the inverse Fourier transform of the calibration region, but these methods are not optimal and need a high number of additional samples. In the following, we will discuss ESPIRiT as a recent subspace-based approach for autocalibration [6].

Calibration can be formulated using the calibration matrix $\mathcal{C}\left\{y_{calib}\right\}$, which is constructed as a block-Hankel matrix from fully sampled calibration data $y_{calib}$ by taking overlapping patches of the multi-channel k-space as rows. Because patches are locally correlated in k-space between receive channels, this matrix has a low rank and singular-value decomposition can be used to compute the basis for the signal subspace $V_{\parallel}$ and for the noise subspace $V_{\perp}$ of patches in k-space.

Given a calibration matrix, highly accurate sensitivities can be recovered using the ESPIRiT algorithm: The null-space condition $V_{\perp}V_{\perp}^{H}R_{r}\hat{y}=0$ for patches $R_{r}\hat{y}$ of the unknown k-space $\hat{y}$ ($r$ enumerating all patches) yields an overdetermined system of equations which can be solved in the least-squares sense using the normal equations:

[TABLE]

This can be rewritten as a a matrix-valued convolution $\mathcal{W}$ that reproduces an ideal fully-sampled k-space ($M$ is the number of samples in a patch):

[TABLE]

Combining this equation $\mathcal{W}y=y$ with the ideal signal model $y=\mathcal{F}S\rho$ and applying an inverse Fourier transform yields $\mathcal{F}^{-1}\mathcal{W}\mathcal{F}S=1S$ where $\rho\neq 0$. In other words, the coil sensitivities appear as eigenvectors to the eigenvalue one of the operator $\mathcal{F}^{-1}\mathcal{W}\mathcal{F}$. As this operator is point-wise in the image domain, it can be efficiently be computed. Interestingly, in the case of certain inconsistencies in the data the algorithm produces multiple sets of eigenvector maps to the eigenvalue one. This can be exploited for more robust parallel MRI by using an extended optimization problem for image reconstruction:

[TABLE]

Here, $S^{i}_{j}$ correspond to multiple sets $c^{i}_{j}$ of sensitivity maps, and $x^{i}$ are multiple sets of images which then have to be combined in a post-processing step.

III Calibrationless Parallel Imaging

For random or non-Cartesian sampling and when a fully sampled calibration region is not available, calibrationless methods are required. For example, exploiting the low rankness of the calibration matrix $\mathcal{C}\left\{y\right\}$ missing k-space samples can be recovered using structured low-rank matrix completion [7]:

[TABLE]

A solution of this problem can be approximated using Cadzow’s algorithm [8], although this is computationally expensive because it requires a large SVD in each iteration. Reconstruction can also be understood directly as a bilinear inverse problem related to multi-channel blind deconvolution [9, 10, 11]. The regularized non-linear inversion (NLINV) algorithm [11] jointly estimates the image $x$ and the coil sensitivities $c_{j}$ based on the formulation as a non-linear inverse problem:

[TABLE]

Here, $W$ a weighting matrix that penalizes high Fourier coefficients to exploit the fact that sensitivities are smooth functions. Fig. 3 shows an experimental example. Inspired by ESPIRiT, this algorithm can be extended to multiple images and sets of sensitivities [12]:

[TABLE]

By introducing new transformed variables $u$ and $v$, and a linear operator $\mathcal{A}$ which maps the outer product $X=uv^{T}$ to $P\mathcal{F}\{c_{j}x\}$, it can be shown that this relaxed version of the algorithm is related to a (convex) structured low-rank matrix recovery problem with nuclear-norm regularization [13].

IV Conclusion

While parallel MRI with known sensitivities is a linear inverse problem for which results from sampling theory and efficient algorithms are readily available, autocalibrating parallel imaging is related to multi-variate multi-channel blind deconvolution and still an active area of research.