## Parallel MRI Reconstruction by Convex Optimization

Cishen Zhang and Ifat-Al Baqee Faculty of Engineering and Industrial Sciences, Swinburne University of Technology Hawthorn, VIC 3122, Australia Email: cishenzhang@swin.edu.au and ibaqee@swin.edu.au

Abstract -In parallel magnetic resonance imaging (pMRI), to find a joint solution for the image and coil sensitivity functions is a nonlinear and nonconvex problem. A class of algorithms reconstruct sensitivity encoded images of the coils first followed by the magnitude only image reconstruction, e.g. GRAPPA. It is shown in this paper that, if only the magnitude image is reconstructed, there exists a convex solution space for the magnitude image and sensitivity encoded images. This solution space enables formulation of a regularized convex optimization problem and leads to a globally optimal and unique solution for the magnitude image reconstruction. Its applications to invivo MRI data sets result in superior reconstruction performance compared with other algorithms.

[6]-[8]. Because of the cross product terms of the image and sensitivity functions, the formulated optimization problem is nonlinear and nonconvex and can only yield local solutions.

Index Terms -Medical imaging; Parallel MRI; Convex optimization; Regularized optimization; Global solution.

## I. INTRODUCTION

Magnetic resonance imaging (MRI) is an advanced modality for obtrusive medical diagnosis which provides very safe scanning, high spatial resolution and flexible contrast for displaying body proton mass. To reduce the duration for scanning without compromising the image quality has been an important and challenging problem in the MRI research. One approach is to implement multiple receiver coils to accelerate MRI scans by acquiring simultaneously undersampled k -space data which is known as parallel MRI (pMRI).

The pMRI reconstruction using undersampled k -space data requires knowledge of spatial sensitivity functions of the multiple receiver coils, which are complex valued and determined by the coil instrumentation the imaged object [1]. Several algorithms have been developed in past years for pMRI reconstruction, one set of algorithms pre-estimate the complex valued sensitivity functions and use the estimated results to reconstruct the the complex valued image such as SENSE [2] and its extensions. Another set of algorithms estimate the sensitivity encoded images of each receiver coil first followed by a image reconstruction operation. Not only knowledge on sensitivity functions is implicitly prerequisite for these algorithms but also the reconstructed image function contains magnitude only values. Typical algorithms of this set are GRAPPA [3] and its extensions using the sum-of-squares (SOS) operation [1]. Some recent algorithms of this group reconstruct the sensitivity encoded images by regularized optimization, e.g. [4], [5]. The third set of algorithms formulate the pMRI reconstruction into regularized optimization problem and jointly estimate the image and sensitivity functions, e.g.

In our previous work [9], it has been shown that, if only the magnitude of the image is reconstructed, the pMRI reconstruction can be formulated into a two-step convex optimization problem which resulted a global optimal solution. To ease the computational burdens of two step formulation, in this paper we present a single step iterative optimization formula for the problem which can lead to efficient computing of the solution as well as yield global optimal solution, also can outperform other pMRI reconstruction algorithms. Without loss of popularity and as the second group of algorithms including GRAPPA and its extensions have done, the magnitude only image reconstruction can meet the needs of most clinical applications.

In this paper, R R , + and C denote the sets of real, nonnegative real and complex numbers, respectively. 〈· , ·〉 denotes the inner product of vectors. glyph[precedesequal] and glyph[followsequal] denote the elementwise operations of ≤ and ≥ on vectors, respectively. glyph[circledot] denotes the Hadamard or elementwise product of vectors. | · | takes elementwise magnitude of vectors and ∠ unitizes elements of vectors, such that v = | v |glyph[circledot] ∠ v for a complex valued vector v . r = ( x, y ) and k = ( k , k x y ) denote the 2D coordinate systems of the spatial image and k -space domains, respectively.

## II. THEORY

## A. Formulation of undersampled k -space data

Consider a pMRI scanner implemented with L receiver coils. Let h ( r ) ∈ C be the 2D spatial MR image function and s i ( r ) ∈ C , i = 1 , · · · , L , be the 2D spatial sensitivity functions of the coils. The sensitivity encoded image functions of the coils are products of h ( r ) and s i ( r ) written as

$$z _ { i } ( { \mathbf r } ) = h ( { \mathbf r } ) s _ { i } ( { \mathbf r } ) \in \mathbb { C }, \ i = 1, \cdots, L.$$

The MRI scan creates the following k -space functions of the L receiver coils

$$\mathfrak { z } _ { \mathfrak { s } } ^ { \mathfrak { t } } \, \quad g _ { i } ( \mathbf k ) = \, \int \, \int z _ { i } ( \mathbf r ) e ^ { - j 2 \pi \langle \mathbf k, \mathbf r \rangle } d \mathbf r, \ \ i = 1, \cdots, L, \ \ ( 1 )$$

which are the Fourier transforms of z i ( r ) = h ( r ) s i ( r ) . The k -space functions are undersampled and formulated into the following discrete data vectors

$$g _ { i } = \mathbf F z _ { i } = \mathbf F ( \mathbf s _ { i } \odot \mathbf h ), \ \ i = 1, 2, \cdots, L, \quad ( 2 )$$

where g i ∈ C M , M &lt; N 2 , i = 1 , · · · , L , are the undersampled k -space data vectors, h , s i , z i ∈ C N 2 are the discretized vectors of h ( r ) , s i ( r ) and z i ( r ) , respectively, and the matrix F ∈ C M N × 2 is the corresponding undersampled partial 2D discrete Fourier transform (DFT) matrix operating on vectors.

For the purpose of magnitude image reconstruction, let h m = | h | ∈ R N 2 + be the magnitude of the image vector h and ˆ s i = s i glyph[circledot] ∠ h ∈ C N 2 + , i = 1 , · · · , L . The undersampled k -space data vectors g i can be written as

$$g _ { i } = \mathbf F z _ { i } = \mathbf F ( \hat { \mathbf s } _ { i } \odot \mathbf h _ { m } ), \ i = 1, 2, \cdots, L. \quad ( 3 )$$

The objective of magnitude image reconstruction is to find a solution for h m using given undersampled k -space data vectors g i , i = 1 , · · · , L .

## B. The convex solution space

Given the undersampled k -space data vectors g i , i = 1 , · · · , L , in the form (2), to find a joint solution for the image vector h and the sensitivity functions s i is in general a nonlinear and nonconvex problem. If only the magnitude of the image function is considered, it is possible to construct a convex solution space for the vector equation (3), which can further lead to a formulation of a convex optimization problem for the magnitude image reconstruction. An observation of the convex solution space is introduced below.

In each bilinear equation z i = ˆ s i glyph[circledot] h m of the sensitivity encoded image functions, there are in general two independent variable vectors which, if known, can determine the third vector variable. Since the sensitivity functions s i have bounded magnitudes due to bounded inductances of the coils, there exist constant vectors b i ∈ R N 2 + such that | ˆ s i | = | s i | glyph[precedesequal] b i , i = 1 , · · · , L . It follows that the magnitudes of z i =ˆ s i glyph[circledot] h m are constrained by

$$| \mathbf z _ { i } | \preceq \mathbf b _ { i } \odot \mathbf h _ { m }, \ i = 1, \cdots, L. \quad \quad ( 4 ) \quad \.$$

If h m and z i are considered as independent variables, for some i , the inequality (4) forms a cone shaped convex hull, which contains the true solutions of h m and z i for the vector equation (3), if the constant bound vector b i is properly chosen. For the scalar case of h m ∈ R + and z i ∈ C , such a convex solution space for h m and z i is displayed in Fig. 1, on top of the complex plane of z i . This convex solution space provides a basis for the convex optimization of the pMRI reconstruction problem and its extension to the high dimensional convex solution space is straightforward. It is, however, noted that the convex solution space only exists for the positive valued magnitude image h m but not for other real or complex valued image vectors.

## C. Convex optimization for image reconstruction

With h m and z i as variables and with given g i and the constraints (4), i = 1 2 , , · · · , L , the k -space data equations (3) are underdetermined and have in general infinite number of solutions for h m and z i . Specifically, if { h m 0 , z i 0 , i = 1 2 , , · · · , L } is a solution for (3), { a h m 0 , z i 0 , i = 1 2 , , · · · , L } is also a solution for any a &gt; 1 . An effective approach

Fig. 1. Convex solution space for h m and z i .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0622v1-figure-c2da0399-3550-44cd-9666-d0a711157a18.jpg)

to finding a meaningful solution for the underdetermined equation (3) is to apply regularized optimization to penalize h m in terms of some normed function of h m .

Let W denotes a general wavelet transformation matrix and assume that the wavelet transformed image vector Wh m is sparse. The glyph[lscript] 1 regularized optimization for the image reconstruction promotes the regularization variable and is formulated as

$$\inf _ { \mathrm n t } \min _ { \mathrm h _ { m }, \mathbf z _ { i }, i = 1, \cdots L } \frac { 1 } { 2 } \sum _ { i = 1 } ^ { L } \| \mathbf g _ { i } - \mathbf F z _ { i } \| _ { 2 } ^ { 2 } + \alpha \| \mathbf W \mathbf h _ { m } \| _ { 1 }$$

$$\text{IRD} \quad \text{subject to} \colon \ h _ { m } \succeq 0, \ \ | z _ { i } | \preceq b _ { i } \odot h _ { m }, \ i = 1, \cdots, L, \\ \dots \quad. \ \. \ \. \ \. \ \. \ \. \ \. \ \. \.$$

where α ≥ 0 is a regularization parameter. The glyph[lscript] 1 regularized optimization problem in (5) can alternatively be formulated as

$$\min _ { \mathbf h _ { m }, z _ { i }, \mathbf w _ { i }, \ i = 1, \cdots L } \ \frac { 1 } { 2 } \sum _ { i = 1 } ^ { L } \| g _ { i } - \mathbf F z _ { i } \| _ { 2 } ^ { 2 }$$

$$+ \alpha \| \text{Wh} _ { m } \| _ { 1 } + \beta \sum _ { i = 1 } ^ { L } \| | z _ { i } | - b _ { i } \odot \text{h} _ { m } + \text{q} _ { i } \| _ { 2 } ^ { 2 } \quad ( 6 )$$

$$\text{subject to} \colon \text{ h} _ { m } \succeq 0, \text{ q} _ { i } \succeq 0, \text{ } i = 1, \cdots, L,$$

where β ≥ 0 is a regularization parameter.

The above formulated regularized optimization problem is convex and can yield a globally optimal and unique solution for the magnitude image h m and the sensitivity encoded image functions z i , i = 1 , · · · , L . It is noted that the glyph[lscript] 1 regularized optimization promotes the sparsity of the regularization variable. It is possible to select other regularization terms of h m and z i to meet reconstruction specifications when needed.

## III. METHODS

## A. k -space data sets

The proposed convex optimization method is applied to two sets of in-vivo MR data for reconstruction and performance evaluation. The first brain data set was acquired on a 3T scanner (GE Healthcare, Waukesha, WI) using an 8 -channel head channel (Invivo, Gainesville, FL) and a 2D T 1 -weighted spin echo protocol with axial plane, TE/TR = 11 700 / ms, FOV= 22 cm 2 , 10 slices and 256 × 256 matrix size. The second spine data set was acquired using an 8 -channel cervicalthoracic-lumbar spine array and a fast spoiled gradient-echo sequence with TR/TE = 300 12 / ms, RBW = 62 5 . kHz, 256 × 256 , Tip angle = 15 ◦ and FOV = 32 × 32 cm 2 . The fully sampled k -space data sets were in the cartesian coordinate system and were manually undersampled by uniform sampling with additional auto-calibration signal (USACS) lines in the phase encoding direction.

## B. Computing Set Ups

The computation of the proposed optimization problem (6) is implemented with the split Bregman method for iterative l 1 regularized optimization in [10]. The regularization parameters are empirically chosen and a tolerance value of 10 -5 is selected for each step of iteration. The haar wavelet matrix is adopted for W . The initial values of h m and z i , i = 1 , · · · , L , are randomly chosen in each iterative reconstruction. And the algorithm is programmed with Matlab (Math-Works, Natick, MA, USA).

To evaluate the reconstruction performance, the reconstructed images, denoted by h o , are compared with the sum of square (SOS) image, which is reconstructed using the fully sampled data and denoted as h SOS . The normalized mean square error (NMSE) of h o is defined as

$$e _ { N M S E } = \frac { \| \mathbf h ^ { o } - \mathbf h _ { S O S } \| _ { 2 } ^ { 2 } } { \| \mathbf h _ { S O S } \| _ { 2 } ^ { 2 } }.$$

The NMSE values of reconstructed images by the proposed algorithm are computed and compared with that by other recently developed methods for the in-vivo brain data sets under the same data reduction conditions.

Fig. 2. Reconstructed brain images by our method, (a) SOS image as reference, (b) Reconstructed image with f nom = 4 , (c) Reconstructed image with f nom = 8 , (d) Reconstructed image with f nom = 12 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0622v1-figure-6e9f2a52-5bfe-4748-9b9b-7f5d25c43684.jpg)

Fig. 3. Comparison of brain image reconstructions with f nom = 8 , (a) by our proposed method, (b) by L 1 CG SPIRiT, (c) by Sparse BLIP, (d) by IRGN TGV.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0622v1-figure-af1e6884-2a2f-41e8-b933-396e09bddc94.jpg)

TABLE I NMSE'S OF DIFFERENT ALGORITHMS

| Method     |   IRGN-TGV |   Sparse-BLIP |   SPIRiT |   Our Method |
|------------|------------|---------------|----------|--------------|
| f nom = 4  |     0.0042 |        0.004  |   0.0066 |       0.0038 |
| f nom = 8  |     0.0056 |        0.0048 |   0.0102 |       0.0045 |
| f nom = 12 |     0.0086 |        0.0072 |   0.0148 |       0.0063 |

## IV. RESULTS

For the brain data set, the reconstructed images by the proposed algorithm using undersampled k -space data of nominal reduction factor, f nom = 4 ( netreductionfactor, f net = 2 66) . and higher f nom = 8 ( f net = 4 00) . and f nom = 12 ( f net = 4 83) . are shown in Fig.2 (b), (c) and (d), respectively, in comparison with the SOS image using the full k space data in Fig.2 (a). For the brain date set of f nom = 8 , Fig.3 presents the reconstructed image by the proposed method in comparison with that by other recent algorithms, which are IRGN-TGV [7], SPIRiT [4] and SparseBLIP [8]. To provide more detailed reconstruction information for comparison, a selected area of each reconstructed image by different algorithms is zoomed and displayed in Fig. 4.

The NMSEs of the reconstructed brain images of different algorithms are listed in Table I for f nom = 4 , 8 , 12 , respectively. The results clearly show that the proposed convex optimization method outperforms other methods in terms of the reconstruction accuracy at high undersampling rates as well as reduction of artifacts.

For the spine data set, the reconstructed image by SOS is given in Fig.5 (a) followed by that reconstructed using the proposed method Fig. 5 (b), for nominal undersampling rate f nom = 4 ( f net = 2 66) . . The NMSE of this reconstructed

Fig. 4. Comparison of a zoomed portion of brain image reconstructions with f nom = 8 , (a) SOS image as reference, (b) Zoomed portion of SOS, (c) by our proposed method, (d) by SPIRiT, (e) by Sparse BLIP and (f) by IRGN-TGV.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0622v1-figure-7cea7d65-b3d4-4229-a164-f2c7359cab2f.jpg)

image is e NMSE = 0 0036 . . For higher reduction rates of f nom = 8 ( f net = 4 00) . and f nom = 12 ( f net = 4 83) . , the corresponding reconstructed images by the proposed algorithm are shown in Fig.5 (c) and 5 (d) with NMSE values 0 0045 . and 0 0061 . , respectively.

With randomly selected different initial values for the magnitude image and sensitivity encoded images for the proposed convex regularization algorithm, the reconstructed brain and spine images with different undersampling rates always converge to the same solutions. These results demonstrated the globally unique and optimal solution of the proposed convex optimization. In contrast, simulations of other nonconvex algorithms resulted in significantly different reconstruction solution under differently initial conditions.

## V. CONCLUSION

The pMRI reconstruction based on undersampled k -space data by optimization methods is a nonlinear and nonconvex problem. The existing optimization methods for jointly computing the image and sensitivity functions can only provide local but not global solutions. The computing of such a nonlinear and nonconvex problem involve complicated procedures and

Fig. 5. Reconstructed spine images by our proposed method, (a) SOS image as reference, (b) Reconstructed image with f nom = 4 , (c) Reconstructed image with f nom = 8 , (d) Reconstructed image with f nom = 12 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0622v1-figure-82d47312-60b7-45ee-b6be-8c03bf6577ea.jpg)

iterations. It is shown in this paper that, if only the magnitude image reconstruction is considered, there is a convex solution space for the magnitude image and sensitivity encoded image functions. This enables formulation of a regularized convex optimization problem for the magnitude image reconstruction and the solution of this problem is globally unique and optimal. Extensive simulations of in-vivo data sets have been carried out and the superior performance of the proposed convex optimization for pMRI reconstruction has been shown.

## REFERENCES

- [1] P. B. Roemer, W. A. Edelstein, C. E. Hayes, S. P. Souza, and O. M. Mueller, 'The NMR Phased Array,' Magnetic Resonance in Medicine , vol. 16, pp. 192-225, 1990.
- [2] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, 'SENSE: Sensitivity encoding for fast MRI,' Magnetic Resonance in Medicine , vol. 42, pp. 952-962, 1999.
- [3] M. Griswold, P. Jakob, R. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, and A. Haase, 'Generalized Autocalibrating Partially Parallel Acquisitions (GRAPPA),' Magnetic Resonance in Medicine , vol. 47, pp. 1202-1210, 2002.
- [4] M. Lustig and J. M. Pauly, 'SPIRiT: Iterative Self-consistent Parallel Imaging Reconstruction From Arbitrary k-Space,' Magnetic Resonance in Medicine , vol. 64, pp. 457-471, 2010.
- [5] S. Park and J. Park, 'Adaptive self-calibrating iterative GRAPPA reconstruction,' Magnetic Resonance in Medicine , vol. 67, pp. 1721-1729, 2012.
- [6] M. Uecker, T. Hohage, K. Block, and J. Frahm, 'Image reconstruction by regularized nonlinear inversion - Joint estimation of coil sensitivities and image content,' Magnetic Resonance in Medicine , vol. 60, pp. 674682, 2008.
- [7] F. Knoll, C. Clason, K. Bredies, M. Uecker, and R. Stollberger, 'Parallel imaging with nonlinear reconstruction using variational penalties,' Magnetic Resonance in Medicine , vol. 67, pp. 34-41, 2012.
- [8] H. She, R. Chen, D. Liang, E. DiBella, and L. Ying, 'Simultaneous image reconstruction and sensitivity estimation in parallel MRI using blind compressed sensing,' 2012, pp. 876 - 879 , 9th IEEE International Symposium on Biomedical Imaging (ISBI) .
- [9] C. Zhang and I. Baqee, 'Parallel magnetic resonance imaging reconstruction by convex optimization ,' 2013, pp. 473-478, 2013 Third International Conference on Innovative Computing Technology (INTECH)

.

- [10] T. Goldstein and S. Osher, 'The split Bregman method for l 1 regularized problems,' SIAM Journal of Imaging Science , vol. 2, pp. 323-343, 2009.
Cishen Zhang Ifat-Al Baqee
Parallel MRI Reconstruction by Convex Optimization
In parallel magnetic resonance imaging (pMRI), to find a joint solution for
the image and coil sensitivity functions is a nonlinear and nonconvex problem.
A class of algorithms reconstruct sensitivity encoded images of the coils first
followed by the magnitude only image reconstruction, e.g. GRAPPA. It is shown
in this paper that, if only the magnitude image is reconstructed, there exists
a convex solution space for the magnitude image and sensitivity encoded images.
This solution space enables formulation of a regularized convex optimization
problem and leads to a globally optimal and unique solution for the magnitude
image reconstruction. Its applications to in-vivo MRI data sets result in
superior reconstruction performance compared with other algorithms.