A Convex Optimization Approach to pMRI Reconstruction
In parallel magnetic resonance imaging (pMRI) reconstruction without using pre-estimation of coil sensitivity functions, one group of algorithms reconstructs sensitivity encoded images of the coils first followed by the magnitude image reconstruction, e.g. GRAPPA. Another group of algorithms jointly computes the image and sensitivity functions by regularized optimization which is a non-convex problem with local only solution. For the magnitude image reconstruction, this paper derives a reconstruction formulation, which is linear in the magnitude image, and an associated convex hull in the solution space of the formulated equation containing the magnitude image. As a result, the magnitude image reconstruction for pMRI is formulated into a two-step convex optimization problem, which produces a globally optimal solution. An algorithm based on split-bregman and nuclear norm regularized optimizations is proposed to implement the two-step convex optimization and its applications to phantom and in-vivo MRI data sets result in superior reconstruction performance compared with existing algorithms.
Magnetic resonance imaging (MRI) is an advanced modality for noninvasive medical diagnosis with continuously growing clinical applications. In comparison with other medical imaging modalities, such as x-ray computed tomography (CT) and ultrasound, MRI is beneficiary in a way that it provides very safe scanning, high spatial resolution and flexible contrast for displaying body tissues. It is, however, also known that the data acquisition of MRI is a relatively long process, which is governed by the time required for physical excitation and relaxation of the radio frequency (RF) magnetic field. The relatively low speed of MRI scan can be an uncomfortable experience for patients and can result in low patient throughput of MRI scan operation. It can also cause motion artifacts in images, typically for respiratory or cardiac organs which have motions during the examination.
To accelerate the scan speed without compromising the image quality has been an important and challenging problem in the MRI research. For this purpose, modern MRI scanners implement multiple receiver coils in phased array mode for parallel acquisition of MRI data in the -space. This technology is known as parallel MRI (pMRI). In pMRI, distinct spatial sensitivities of the receiver coils can enable simultaneous acquisition of -space data containing complementary information. As a result, the pMRI can accelerate MRI scans to considerably reduce the overall scan time and a combined set of undersampled partial -space data acquired from different receiver coils can provide sufficient information for image reconstruction.
The MRI produces gray value images displaying the spatial magnetic spin density function of the imaged object. The magnetized spin density function is complex valued with magnitude and phase and is determined by the proton density of the imaged object, the external magnetic field and RF excitation pulses. The coil sensitivity functions are also complex valued and bounded due to their finite inductance values . To reconstruct the MR image from the -space data of a single receiver coil scanner is an inverse Fourier transform process under the uniform sensitivity assumption. But the pMRI reconstruction using undersampled -space data is not a straightforward task. It requires knowledge of spatial sensitivity functions of the multiple receiver coils, which are in general not only determined by the coil instrumentation but also dependent upon the imaged object. There have been numerous algorithms developed in past years for pMRI reconstruction. Depending on how the information of sensitivity functions is incorporated into the image reconstruction and how the image information is reconstructed, the existing reconstruction algorithms can be classified into three groups.
The first group of algorithms pre-estimates the complex valued sensitivity functions and use the estimated results to reconstruct the magnitude and phase functions of the complex valued image. The performance of the algorithms depends on the accuracy of the pre-estimated sensitivity functions. Typical algorithms of this group are SMASH, SENSE , and their extensions such as [4, 5, 6]. Also included in this group are some algorithms based on estimation of sensitivity functions and regularized optimization e.g. [7, 8, 9, 10, 11]. The second group of algorithms estimates the sensitivity encoded images of each receiver coil first followed by a image reconstruction operation to obtain the image. These algorithms do not require knowledge of the sensitivity functions but the reconstructed image function is magnitude only without containing the phase information. Typical algorithms of this group is GRAPPA , IIR GRAPPA  and their extensions using the sum-of-squares (SOS) operation . There are also recent algorithms which reconstruct the sensitivity encoded images by regularized optimization, e.g. [14, 15, 16, 17]. The third group of algorithms formulates the pMRI reconstruction into a regularized optimization problem without requiring estimation of sensitivity functions [18, 19, 20, 21, 22]. The algorithms jointly compute the complex valued image and sensitivity functions by minimizing a performance index function which incorporates the reconstruction error and regularization terms. Because of the inherent cross product terms of the image and sensitivity functions, the formulated optimization problem is bilinear in the optimization variables and hence nonconvex. It therefore can only result in a local solution depending on the selection of initial condition and has difficulties in finding the global optimal solution as well as computational complexity.
This paper considers the pMRI reconstruction problem by regularized optimization without using knowledge of the sensitivity functions and tackles the difficulties of nonconvex optimization and local solution of the third group of algorithms. It is shown that, if only the magnitudes of the image and sensitivity functions are reconstructed, the pMRI reconstruction can be formulated into a linear and convex optimization problem which has a global optimal solution. The linear and convex formulation of the problem can lead to efficient computing of the solution and the global optimal solution 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. It is also noted that there are some application cases, such as phase contrast imaging for detection of the velocity of flow , where the phase infromation of the image is required and the magnitude only image reconstruction is not sufficient.
Like the two-step procedure of GRAPPA and its extensions, which first estimate the sensitivity encoded image functions of each coil followed by an SOS operation to construct the magnitude image, the proposed pMRI reconstruction in this paper is formulated into a two-step convex optimization problem, with the first step optimization solving the sensitivity encoded image functions of each coil and the second step optimization solving the magnitude image function. The two optimization steps operate sequentially in a way that the second step optimization is carried out after completion of the first step optimization, which is different from the iterative alternating optimization for solving the nonlinear optimization problems. The two-step convex optimization is implemented with an algorithm based on Split-bregman method and nuclear norm regularization and applied to in-vivo -space data for pMRI reconstruction. Its performance in reconstruction accuracy and efficiency in comparison with other methods is demonstrated.
In this paper, , and denote the sets of real, nonnegative real and complex numbers, respectively. The lower bold case letter denotes vectors and the capital bold case letter denotes matrices. and denote the elementwise operations of and on vectors, respectively. denotes the Hadamard or elementwise product of vectors. takes elementwise magnitude of vectors and unitizes elements of vectors, such that for a complex valued vector . denotes the inner product of vectors. and denote the 2D coordinate systems of the spatial image and -space domains, respectively.
Ii Formulation of the -space data for convex optimization
Ii-a The undersampled space data
Consider a pMRI scanner implemented with receiver coils. Let be the 2D spatial MR image function and , , be the 2D spatial sensitivity functions of the coils. The sensitivity encoded image functions of the coils are , , which are products of and . The MRI scan creates the following -space functions of the receiver coils
which are the Fourier transforms of .
The discrete version of the Fourier transform equations in (1) can be written in the vector forms as
where , , are the discretized vectors of , , and , respectively, and the matrix operates the 2D discrete Fourier transform (DFT) on the vectorized form of 2D matrices. The undersampled vectors of the -space data , denoted by , with , can be represented as
where is the corresponding undersampled 2D DFT matrix operating on vectors.
Ii-B The convex solution space
Given the undersampled -space data vectors in the form (3), to find a joint solution for the image and sensitivity functions 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 magnitude image and the sensitivity encoded functions , which can further lead to a convex optimization formulation of the image reconstruction. An intuitive observation of the convex solution space is provided below.
Let be the magnitude of the image vector . Since the magnitudes of are bounded due to bounded inductances of the coils, there exist constant vectors such that , . It follows that
In each bilinear equation of the sensitivity encoded image functions, there are two independent variable vectors which, if known, can determine the third vector variable. If and are considered as the solution variables, the inequalities (4) form a cone shaped convex hull containing the solutions of and , with properly chosen constant bound vectors . Such a convex solution space is displayed in Fig. 1, on top of the complex plane of , for the scalar case of and . This 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 magnitude image but not for any other real or complex valued image vectors.
Ii-C Linear formulation of the pMRI reconstruction
The solution space for and , as shown in Fig.1, displays the convex nature of the problem. But, as seen from (3), the magnitude image is a bilinear variable, coupled with , of the composite vectors , , so is not a linear variable of the problem equation. To facilitate formulation of a linear model for the convex optimization, introduce the magnitude and phase vectors of denoted by and , respectively, such that , . Further introduce decoupling parameter vectors , , and denote their corresponding diagonal matrices as . Using , the magnitude vectors can be written as
, where . It follows from that are linearly and elementwisely bounded by , i.e.
Let , and be the stacked vectors of and , , respectively. The undersampled -space vectors in (3) can be rewritten as
where is the blocked diagonal matrix of . The stacked vector form of in (5) is
The vector equation (8) shows that the magnitude of the image function is linearly decoupled from the magnitude vector of the sensitivity encoded image functions. This technical result is instrumental for the proposed convex optimization for pMRI reconstruction.
Iii Two-step convex optimization for pMRI reconstruction
Iii-a General formulation
It is observed from (7) that the global solution for the sensitivity encoded image vector and hence its magnitude and phase can be obtained by solving the linear equation . Using the solution for , the linear equation (8) can be further solved to obtain a solution for the magnitude of the image. Once the magnitude is obtained. The magnitudes , , of the sensitivity functions can be further determined using (5). Based on this observation, the general formulation of the proposed pMRI reconstruction consists of two sequential convex optimization steps and as follows.
Based on , the first step solves the complex valued sensitivity encoded image vector and hence its magnitude and phase by the following regularized convex optimization
where is a convex regularization function to be further specified according to application conditions.
Suppose that is the optimal solution of the optimization problem with and being the corresponding optimal solutions of and , respectively. Substituting into (8) yields
Given , the second step of the proposed convex optimization is
where is a convex regularization function to be further specified according to application conditions, are the -fold stacked vector of and the stacked vector of , , respectively. Since the linear equation (10) is underdetermined and has infinite number of solutions, the inequality based on (6) forms a convex hull to constrain and in the solution space. The solution for provides the optimal magnitude image as well as the optimal . The corresponding optimal solutions , , for the magnitude of sensitivity functions can be further determined by . At some points where the image function has zero values, feasible solution values for sensitivity functions are not available. Proper interpolation may be introduced at these points for reconstruction of the magnitude sensitivity functions based on their smooth property.
In the two optimization steps and , is to optimize the solution for (7) and is originally a linear problem. The convexity of the optimization problem is built upon the decoupled linear equation (10) and convex solution space specified by (6). This is possible only if the solution variable is the magnitude only image vector. For complex valued and , although a decoupling linear equation in the same form of (10) can be formulated, a convex set in the solution space does not exist for the decoupled variables, so the problem remains nonlinear and nonconvex.
Iii-B Split-bregman and nuclear norm regularized optimizations
The above formulated convex optimization steps and are in general forms and can be implemented with different regularity functions and variable constraints. Taking into account properties of MR images and sensitivity functions, this subsection presents an implementation of with split-bregman method and the nuclear norm regularized optimization for implementation and computation of .
The sensitivity encoded function to be optimized in problem is a product of the image and sensitivity functions and typically can have a piecewise smooth characteristic. The implementation of problem takes this characteristic into account and incorporates it into the regularization function . It is known that Bregman iteration  can solve a broad class of regularized optimization problems. It can result in superior image reconstruction performance when a hybrid of Bounded-Variation (BV) and Besov regularization is used  and has been applied to MR image reconstruction, e.g. . Applying the split-bregman method to , the regularization function can be formulated as
where denotes the bounded variation norm defined as with and being the difference operators in the and directions, respectively, is a wavelet transform matrix and and are regularization parameters. The regularization terms and in (12) are to promote, respectively, the piecewise smoothness of and energy compactness of . The two regularization terms, together with (9), yield the following split-bregman regularized optimization for solving
In the optimization of problem , a reduction of the magnitude in the solution of the underdetermined linear equation (10) can result in the value of and hence the values of to grow. Thus it requires a proper scale of the solutions for and by the regularization function and appropriate constraints on the solutions. The implementation of considers the nuclear norm regularized optimization  which has shown promising results in computational efficiency and accuracy in the application to MR image reconstruction , . The nuclear norm of the magnitude image vector, denoted by , is defined as the sum of singular values of . To use as a regularization term together with the inequalities , which specify that the solution for is linearly bounded by in a convex hull, can provide proper scaling and effective constraints on and for computing the solution for . As a result, the convex optimization problem is formulated as
The above convex optimization problem (14) can be equivalently formulated as
where and are regularization parameters.
Iv Phantom and in-vivo data experiments
Iv-a Cartesian and Non-Cartesian Data
Two sets of in-vivo MRI data were adopted to test the proposed convex reconstruction algorithm. The first is a brain data set (available at http://black.bme.ntu.edu.tw/tool_sense.html) of a healthy human volunteer was acquired by a 3 Tesla SIEMENS Trio scanner with an eight-channel head array and an MPRAGE (3D Flash with IR prep.) sequence. The parameters of the scan were ms, , , flip angle , slice Thickness mm and mm. The fully acquired -space data in the cartesian co-ordinate system are manually undersampled to obtain the uniform sampling with additional auto-calibration signal (USACS) patterns of -, -, - and -fold acceleration rates, respectively, with additional 36 extra auto-calibrating signal (ACS) lines in the central space region along the phase encoding direction in each pattern. As a result, the corresponding net reduction factors are , and , respectively.
Another data set of spine (available at http://ece.tamu.edu/jimji/pulsarweb/downloads.htm) was acquired from a -channel cervical-thoracic-lumbar spine array using a fast spoiled gradient-echo sequence and parameters ms, , , tip angle and cm. The fully acquired -space data in the cartesian coordinate system are undersampled to generate the USACS patterns of -, - and -fold acceleration rates, respectively, with each pattern having extra ACS lines along the phase encoding direction.
The proposed two-step convex optimization algorithm is also applied to a set of scanned non-cartesian phantom data, which is available at http://www.eecs.berkeley.edu/m̃lustig. The phantom data set was scanned on a GE Signa-Excite 1.5T scanner using a -channel cardiac coil set with a spiral gradient echo sequence. The spiral trajectory was designed with interleaves, cm field of view, mm in-plane resolution and readout time of ms. The -space Data was undersampled by choosing out of interleaves. For image reconstruction in case of non-cartesian data sets, the NUFFT code by  was applied.
Iv-B Computational set ups
The computation of the split-bregman based optimization problem (13) was implemented using the iterative algorithm proposed in . The algorithm proposed in , together with the result in , was adopted for resolving the nuclear norm regularized optimization problem (14). LSQR  tools were used in nuclear norm regularized optmization. For wavelet transformations in the Wavelet regularized reconstruction, David Donoho’s Wavlab codes were used. Two wavelet familes, ”Haar” and ”Daubechies” were selected as the sparsifying transform basis. The regularization parameters were empirically chosen and a tolerance value of is selected for each step of iteration. Both algorithms are programmed with Matlab (Math-Works, Natick, MA, USA).
To evaluate the reconstruction accuracy, the reconstructed images, denoted by , are compared with the sum of square (SOS) image, denoted as , which is reconstructed using the fully sampled -space data. The the normalized mean square error (NMSE) of is defined as
The reconstructed images by the proposed algorithm are computed and compared with the reconstructed images by conjugate gradient (CG)-SPIRiT  with penalty, GRAPPA , JSENSE and IRGN-TGV  algorithms for the in-vivo channel brain data sets under the same data reduction conditions. The Matlab codes as well as the regularization parameters and initial conditions, where applicable, for computations of these algorithms are originated from http://www.eecs.berkeley.edu/~mlustig/Software.html for GRAPPA and CG-SPIRiT, http://cai2r.net/sites/default/files/software/irgntv.zip for IRGN-TGV and https://pantherfile.uwm.edu/leiying/www/index_files/software.htm for JSENSE.
The global solutions of the proposed convex optimization algorithm are tested with different initial image conditions. Two typical initial image conditions, a randomly generated matrix and diagonal lines matrix, shown in Fig. 2, of compatible dimensions with reconstructed images.
Iv-C In-vivo cartesian reconstruction
Fig.3 (a), (c), (e) and (g) display the reconstructed images from manually under-sampled data of acceleration factors and , respectively, for the eight-channel brain image in comparison with the reference image reconstructed by SOS from the full -space data set. The regularization parameters are selected as , , and , using the ”Haar” transform. At a smaller acceleration rate such as , the reconstructed image portrays good quality with very small difference from the reference image. Some quality degradation can be observed from the reconstructed images with higher acceleration factor such as . The NMSEs of the images are computed as , , and for the acceleration factors of and , and the corresponding error images are shown in Fig.3 (b), (d), (f) and (h), respectively.
For the brain date set of acceleration factor , Fig.4 presents the reconstruction image of the proposed method in comparison with the images reconstructed by other commonly known algorithms, which are IRGN-TGV, CG-SPIRiT, GRAPPA and JSENSE. A selected area, as blocked within the marked rectangle in the reference image of Fig.4, is zoomed for each reconstructed image and jointly displayed with the corresponding full size image for ease of comparison. Among these algorithms, GRAPPA and CG-SPIRiT are members of the second group using the SOS operation and IRGN-TGV and JSENSE are nonlinear optimization algorithms of the third group. Noticeable errors and artifacts are shown in reconstruction results of Fig.4 (c), (d) and (e) of CG-SPIRiT, GRAPPA and JSENSE algorithms, respectively. More careful observation can also find artifacts in the zoomed image Fig.4 (g) reconstructed by the nonlinear iterative algorithm IRGN-TGV. The NMSE values of the reconstructed images by the different algorithms are listed in Table I. The the average computational time durations of repeatedly running the original Matlab codes of the different reconstruction algorithms on a workstation with Intel Xeon Processor E5-2609 and 16 GB RAM are given in Table II.
For the in-vivo spine data set, the regularization parameters of the proposed algorithm are chosen as , , and . The wavelet transform matrix is ”Haar”. The reference image constructed by the full data set and SOS operation is given in Fig.5, followed that the reconstructed images by the proposed method, for nominal undersampling rate in Fig.5(a). The estimated NMSE of this reconstructed image is . For higher reduction rates of and , the corresponding reconstructed images by the proposed algorithm are shown in Fig.5(b) and 5(c), with estimated NMSE values of and , respectively.
Iv-D Non-cartesian data reconstruction of phantom
For the phantom data of the spiral pattern, a reference image as given in Fig.6 is produced by applying the NUFFT and SOS operations on the full data set. The proposed algorithm, with regularization parameters , , , and the wavelet transform matrix ”Daubechies”, is applied to the undersampled spiral phantom data and produces the reconstructed image in Fig.6(a) and the corresponding error image, with respect to the reference image, in Fig.6(c). Another algorithm CG-SPIRiT which is capable of non-cartesian reconstruction is also applied to the same undersampled data set, resulting in the reconstructed image in Fig.6(b) and the corresponding error image in Fig.6(d). A selected area from both reconstruction results are cropped and scaled up, as shown in Fig.6 (e) and (f) respectively, to demonstrate the effectiveness of our proposed algorithm.
The proposed convex optimization approach to pMRI reconstruction is build upon a convex solution space which exists only for the magnitude image function but does not exist for any real and complex valued images. This paper formulated a two-step convex optimization to solve the pMRI reconstruction and it is possible to solve the convex optimization problem with alternative formulations.
The solution space of the proposed two-step optimization is a convex hull specified by the constraints in (14) which contains the true solution for the image and sensitivity functions. In general, the optimal solution and its computation depend on selections of the regularization parameters as well as the constraint vector in (14). A priori knowledge of the image and sensitivity functions and empirical tests of the parameter and constraints can be helpful for efficient and accurate computing of the solution. The proposed algorithm produces a global solution in the sense that the solution is unique and independent of the initial image value for the computational algorithm. This is a distinctive characteristic of the proposed method, because all other existing methods for optimization of the image reconstruction, without using previously estimated sensitivity values or the SOS operation, can provide only local solutions, which are dependent on the initial value of the algorithm. In the phantom and in-vivo data reconstructions by the proposed method, all global solutions of the proposed algorithms were tested with different initial conditions including that shown in Fig.2 and their uniqueness was verified. In contrast, the solutions of other algorithms based on non-convex optimization, such as IRGN-TGV and JSENSE, are local only. Our experiments showed that their reconstruction results are very different subject to different initial conditions.
The experiments of the in-vivo and phantom data sets demonstrated better image reconstruction quality of the proposed method than that of GRAPPA and CG-SPIRiT which are algorithms of the second group using SOS operation. It indicates that the optimization with properly specified regularization terms can provide better reconstruction results than that of the simple SOS operation. This reconstruction improvement, however, involves more workload in the iterative computing of the optimal solution, which can be seen from TableII of the computational time durations. Because of the linear and convex nature of our proposed algorithm, it has faster and more efficient computation of the optimal solutions in comparison with the nonlinear optimization algorithms IRGN-TGV and JSENSE as shown by their computational time durations in TableII.
The proposed algorithm for computing the phantom and in-vivo images is only a specific realization of the general linear and convex optimization method in terms of the two-step optimization problems and . Algorithms using other regularization terms and realizations relevant to different reconstruction requirements can be possible and will be subject to future studies.
The reconstruction of MR images based on undersampled -space data by optimization methods for pMRI is known as a nonlinear and nonconvex problem. It is a recently active research area in MRI reconstruction and the existing optimization methods without using estimated sensitivity functions or the SOS operation can only provide local but not global solutions. And the solutions for such a nonlinear and nonconvex problem involve complicated computational procedures and iterations. In this paper, a linear equation is derived for the undersampled -space data set in terms of the magnitude image function, which enables the formulation of the pMRI reconstruction into a two-step convex optimization problem. It is applicable to both cartesian and non-cartesian data sets and can provide a globally optimal solution and faster computation for the pMRI reconstruction problem. An algorithm is presented in this paper and applied to phantom and in-vivo MRI data sets to demonstrate the reconstruction performance and effectiveness of the proposed method.
-  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.
-  D. K. Sodickson and W. J. Manning, “Simultaneous acquisition of spatial harmonics (SMASH): Fast imaging with radiofrequency coil arrays,” Magnetic Resonance in Medicine, vol. 38, p. 591â603, 1997.
-  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.
-  W. Kyriakos, L. Panych, D. Kacher, C. Westin, S. Bao, R. Mulkern, and F. Jolesz, “Simultaneous acquisition of spatial harmonics (SMASH): Fast imaging with radiofrequency coil arrays,” Magnetic Resonance in Medicine, vol. 44, p. 301â308, 2000.
-  C. Liu, R. Bammer, and M. Moseley, “Parallel imaging reconstruction for arbitrary trajectories using k-space sparse matrices (kSPA),” Magnetic Resonance in Medicine, vol. 58, p. 1171â1181, 2007.
-  B. Madore, “UNFOLD-SENSE: a parallel MRI method with self-calibration and artifact suppression,” Magnetic Resonance in Medicine, vol. 52, pp. 310–320, 2004.
-  A. Majumdar and R. K. Ward, “Nuclear norm-regularized SENSE reconstruction,” Magnetic Reson Imaging, vol. 30, pp. 213–221, 2012.
-  Y. Chen, W. W. Hager, F. Huang, D. T. Phan, X. Ye, and W. Yin, “Fast Algorithms for Image Reconstruction with Application to Partially Parallel MR Imaging,” SIAM Journal on Imaging Sciences, vol. 5(1), pp. 90–118, 2012.
-  B. Liu, F. M. Sebert, Y. Zou, and L. Ying, “SparseSENSE: randomly-sampled parallel imaging using compressed sensing,” 2008, p. 3154, Proceedings of the 16th Annual Meeting of ISMRM, Toronto.
-  D. Liang, B. Liu, J. J. Wang, and L. Ying, “Accelerating SENSE Using Compressed Sensing,” Magnetic Resonance in Medicine, vol. 62(6), pp. 1574–1584, 2009.
-  F. Huang, Y. Chen, W. Yin, W. Lin, X. Ye, W. Guo, and A. Reykowski, “A Rapid and Robust Method for Sensitivity Encoding with Sparsity Constraints: Self-feeding Sparse SENSE,” Magnetic Resonance in Medicine, vol. 64(4), pp. 1078–1088, 2010.
-  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.
-  Z. Chen, J. Zhang, R. Yang, P. Kellman, L. A. Johnston, and G. F. Egan, “IIR GRAPPA for Parallel MR Image Reconstruction,” Magnetic Resonance in Medicine, vol. 63, pp. 502–509, 2010.
-  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.
-  M. Murphy, M. Alley, J. Demmel, K. Keutzer, , S. Vasanawala, and M. Lustig, “Fast -SPIRiT Compressed Sensing Parallel Imaging MRI: Scalable Parallel Implementation and Clinically Feasible Runtime,” IEEE Transactions on Medical Imaging, vol. 31, pp. 1250–1262, 2008.
-  D. S. Weller, J. R. Polimeni, L. Grady, L. L. Wald, E. Adalsteinsson, and V. K. Goyal, “Sparsity-Promoting Calibration for GRAPPA Accelerated Parallel MRI Reconstruction,” IEEE Transactions on Medical Imaging, vol. 32, pp. 1325–1335, 2013.
-  S. Park and J. Park, “Adaptive self-calibrating iterative GRAPPA reconstruction,” Magnetic Resonance in Medicine, vol. 67, pp. 1721–1729, 2012.
-  L. Ying and J. Sheng, “Joint image reconstruction and sensitivity estimation in SENSE (JSENSE),” Magnetic Resonance in Medicine, vol. 57, pp. 1196–1202, 2007.
-  M. Uecker, T. Hohage, K. T. 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. 674–682, 2008.
-  S. H, C. RR, D. Liang, C. Y, and Y. L, “Image Reconstruction From Phased-Array MRI Data Based on Multichannel Blind Deconvolution,” 2010, pp. 760–763, 7th IEEE International Symposium on Biomedical Imaging: From Nano to Macro.
-  D. Gol and L. C. Potter, “Ambiguity and Regularization in Parallel MRI,” 2011, pp. 2829–2832, 33rd Annual International Conference of the IEEE Engineering-in-Medicine-and-Biology-Society (EMBS).
-  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.
-  N. J. Pelc, R. J. Herfkens, A. Shimakawa, and D. R. Enzmann, “Phase contrast cine magnetic resonance imaging,” Magnetic Resonance Quarterly, vol. 4, pp. 229–254, 1991.
-  T. Goldstein and S. Osher, “The split bregman method for regularized problems,” SIAM Journal of Imaging Science, vol. 2, pp. 323–343, 2009.
-  A. Majumdar and R. K. Ward, “An algorithm for sparse MRI reconstruction by Schatten p-norm minimization,” Magnetic Reson Imaging, vol. 29, pp. 408–417, 2011.
-  L. Bregman, “The relaxation method of finding the common points of convex sets and its application to the solution of problems in convex optimization,” USSR Computational Mathematics and Mathematical Physics, vol. 7, pp. 200–217, 1967.
-  B. Liu, “Regularized Sensitivity Encoding (SENSE) Reconstruction Using Bregman Iterations,” Magnetic Resonance in Medicine, vol. 61, pp. 145–152, 2009.
-  B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed minimum rank solutions to linear matrix equations via nuclear norm minimization,” SIAM Review, vol. 52, pp. 471–501, 2010.
-  J. A. Fessler and B. P. Sutton, “Nonuniform fast Fourier transforms using minmax interpolation,” IEEE Transactions on Signal Process, vol. 51, pp. 560–574, 2003.
-  A. Dax, “A hybrid algorithm for solving linear inequalities in a least squares sense,” Numerical Algorithms), vol. 50, pp. 97–114, 2009.
-  J. Buckheit and D. L. Donoho, Wavelets and Statistics. Springer-Verlag, Berlin, 1996.