Adaptive Transform Learning and Joint Sparsity Based PLORAKS Parallel Magnetic Resonance Image Reconstruction

Parallel magnetic resonance (MR) imaging is an important acceleration technique based on the spatial sensitivities of array receivers. The recently proposed Parallel low-rank modeling of local k-space neighborhoods (PLORAKS) approach uses the low-rank matrix model based on local neighborhoods of undersampled multichannel k-space data for reconstruction purposes. The joint total variation (JTV) regularization term was then combined with the PLORAKS model to improve the quality of reconstructed images. To further improve the quality of parallel MR imaging, we propose combining adaptive transform learning and joint sparsity with the PLORAKS model to obtain two algorithms, and reconstruction problems are solved by using the alternating direction method of multipliers (ADMM) and conjugate gradient techniques. The experimental results show that the two proposed algorithms can achieve higher performance than the PLORAKS algorithm and the PLORAKS-JTV algorithm with the JTV regularization term in terms of the signal-to-noise ratio (SNR), normalized root mean square error (NRMSE), high-frequency error norm (HFEN), and structural similarity index measure (SSIM).


I. INTRODUCTION
Magnetic resonance (MR) imaging is an important imaging tool for modern medical diagnoses because it has no ionizing radiation and can provide good contrast between different soft tissues. However, the scanning speed of MR imaging is fundamentally limited due to physical and physiological constraints. Increasing the scanning speed of MR imaging has been a research hotspot in recent years. Reducing the amount of acquired data is a commonly used method to accelerate MR imaging. Therefore, it is very important to reconstruct a high-quality MR image from the undersampled data. According to the emerging compressed sensing (CS) theory [1], [2], if the signal is sparse in some transform domains, it can likely The associate editor coordinating the review of this manuscript and approving it for publication was Hengyong Yu . be accurately reconstructed the signal from highly undersampled observation data by using nonlinear optimization methods. Since MR images are sparse in some transform domains, such as wavelet transforms and spatial finite differences, they can be reconstructed from highly undersampled k-space data according to CS theory [3]. Some CS-based MR image reconstruction methods have been proposed. The authors in [3]- [6] proposed using the composite regularization terms of 1 norms of wavelet coefficients and total variation (TV) to improve the reconstruction quality of the MR images. The authors in [7], [8] proposed using the patch-based directional wavelets (PBDW, PBDWS) to improve the quality of reconstructed MR images. The authors in [9] proposed integrating a novel iterative feature refinement (IFR) module with CS to restore meaningful structures and details without introducing too much additional complexity. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ In addition to the above methods of fixed sparse representation, researchers have also developed reconstruction methods based on adaptive sparse representation; these include MR image reconstruction algorithms based on dictionary learning [10] and reconstruction algorithms based on sparse transform learning [11]. In dictionary learning based MR image (DLMRI) reconstruction algorithm [10], dictionary learning and reconstruction are performed simultaneously, and k-means-based singular value decomposition (K-SVD) and orthogonal matching pursuit (OMP) are used for dictionary learning and sparse coding, respectively. However, dictionary learning and sparse coding require iteration, and this results in high computational complexity. Therefore, sparse transform learning based reconstruction algorithms for MR images (TLMRI) [11]- [16] have been proposed. In these algorithms, since transform learning and sparse coding have analytical solutions, the computational complexity is reduced.
Parallel MR imaging technology is another method for accelerating MR imaging, that is often combined with compressed sensing technology to improve the quality of reconstructed images. Parallel MR imaging uses multiple coils with different spatial sensitivities to receive data [17], [18]. Some methods, such as simultaneous acquisition of spatial harmonics (SMASH) [19], and sensitivity encoding (SENSE) [20], use the sensitivity information explicitly. The authors in [21]- [24] combined the SENSE model with the TV regularization term, and improved the reconstruction quality. However, it is difficult to measure sensitivity information with high accuracy in practical applications. Small errors in sensitivity can cause artifacts to persist in a reconstructed image. Uecker et al. proposed an ESPIRiT model [25] that could estimate multiple sets of coil sensitivities. Duan et al. [26] proposed combining the ESPIRiT model with the joint TV regularization term, and effectively improved the reconstruction quality. Another class of methods utilizes the sensitivity information implicitly to avoid the difficulty of measuring or estimating sensitivity; these methods include generalized auto-calibrating partially parallel acquisitions (GRAPPA) [27], which requires only a small amount of fully sampled auto-calibration signals (ACSs) for calibration. Iterative selfconsistent parallel imaging reconstruction (SPIRiT) [28] was a GRAPPA-based auto-calibrating method, that enforced calibration consistency for every point (acquired and unacquired) on the grid. Murphy et al. proposed combining the SPIRiT model with the regularization terms of the joint 1 norms of wavelet coefficients, and solving the reconstruction problem by using the projection over convex sets (POCS) algorithm [29]. The authors in [30], [31] proposed combining the SPIRiT model with the joint TV regularization term, and solving reconstruction problems by using the alternating direction method of multipliers (ADMM) and the operator splitting (OS) technique, respectively.
Recently, several calibrationless parallel imaging reconstruction methods have been proposed. The calibration-free locally low-rank encouraging reconstruction (CLEAR) method [32] uses the locally low rank (LLR) matrix constraint of the image domain for reconstruction, and Saucedo et al. proposed improving the computational efficiency of LLR based parallel MR imaging reconstruction [33]. The authors in [34], [35] proposed using the joint sparse model of the image domain for the reconstruction process. The joint TV regularization term was used for the reconstruction process in [36]. Both simultaneous autocalibrating and k-space estimation (SAKE) [37] and PLORAKS methods [38] use the low-rank matrix model based on the local neighborhoods of undersampled multichannel k-space data. Haldar et al. also proposed combining the PLORAKS model with the joint TV regularization term to further improve the reconstruction quality [38]. To further improve the reconstruction accuracy of parallel MR imaging, the authors in [39] proposed learning joint-sparse codes for calibrationfree parallel MR imaging (LINDBERG) by formulating the parallel MR imaging problem as an 2 − F − 2,1 minimization objective.
In addition, researchers have proposed algorithms based on deep learning for parallel MR imaging and dynamic MR imaging. Specifically, Wang et al. proposed a dynamic MR imaging method with both k-space and spatial prior knowledge integrated via multi-supervised network training (DIMENSION) [40]; the authors in [41] exploited a residual complex convolutional neural network for fast parallel MR imaging.
In this paper, we propose two new parallel MR image reconstruction methods to further improve the quality of reconstructed parallel MR images by combining joint sparsity and adaptive transform learning with the PLORAKS model. In summary, our contributions can be summarized as follows: Since adaptive patch-based transforms can capture local image features effectively, we propose combining adaptive transform learning with the PLORAKS model to improve the quality of reconstructed parallel MR images. To make full use of the correlations between multicoil images, we propose to introduce the joint sparsity of multiple coils. And then, reconstruction problems are solved by using the ADMM and conjugate gradient (CG) techniques. Experimental results show that the proposed methods are superior to PLORAKS and PLORAKS-JTV, and these results demonstrate the effectiveness of this method in parallel MR imaging.
The structure of this paper is as follows: In Section II, we review the LORAKS and PLORAKS models. Then, we combine joint sparsity and adaptive transform learning with the PLORAKS model, and we present our algorithms (PLORAKS-TL and PLORAKS-JTL) for solving optimization problems in Section III. Simulated experimental results and analyses are provided in Section IV. Finally, Section V concludes this paper.

II. OVERVIEW OF LORAKS AND PLORAKS
The LORAKS model was a new MR image reconstruction model based on the fact that low rank matrices could be constructed from fully-sampled k-space data when the given 212316 VOLUME 8, 2020 image had limited support or a smooth phase. The reference work [42] described three matrices: the C matrix generated from the linear dependence of the image, the S matrix, and the G matrix constructed using the smoothly varying spatial image phase. We employ the S matrix for reconstruction due to its high reconstruction quality.
The S matrix [38], [42] can be constructed by a linear operator . The definition of R S (·) is as follows: where S r + , S r − , S i + , and S i − ∈ R K ×N R are expressed as: The undersampled k-space data can be expressed as: where f ∈ C Q is the fully-sampled k-space data, d ∈ C M is the undersampled k-space data, and A ∈ C M ×Q is the undersampling matrix. The LORAKS reconstruction problem can be formulated as the following optimization problem: where λ represents the regularization parameter, and J r (X ) is a regularization function. In reference work [42], J r (X ) is defined by: where X ∈ R 2K ×2N R , and σ k is the k th singular value of matrix X . PLORAKS extends the LORAKS model to parallel MR image reconstruction. In parallel MR imaging, f 1 , . . . , f l , . . . , f L ∈ C Q are fully-sampled k-space data from L different receiver coils, and f P = [f 1 , . . . , f l , . . . , f L ] ∈ C LQ . In PLORAKS, the modified S matrix can be constructed by a linear operator R PS (·) : C LQ → R 2K ×2LN R , which is defined by: The PLORAKS reconstruction problem can be expressed as the following optimization problem: where d P = [d 1 , . . . , d l , . . . , d L ] ∈ C LM , and d l ∈ C M is the undersampled k-space data of the l th coil.
Reference [38] also combined the joint TV (JTV) regularization term with the PLORAKS model, and the corresponding reconstruction problem is as follows: where F represents the coil-by-coil two-dimensional Fourier transform, · JTV represents the JTV regularization term, and α is the regularization parameter. For convenience of expression, we omit the superscript 'P' for the variables in the following description.

III. THE PROPOSED ALGORITHMS A. PROBLEMS FORMULATION
Although the PLORAKS model combined with the JTV regularization term can effectively improve the quality of reconstructed images, there is still much room for improvement. To further improve the reconstruction quality of parallel MR imaging, this paper combines adaptive transform learning [11], [13] with the PLORAKS model and obtains the following optimization problem: where P j (·) : C LQ → R n×L is the j th linear operator that extracts an image patch of size √ n× √ n from each coil image, and then converts these L image patches into a C n×L matrix. W ∈ R n×n represents an adaptive sparse transform learned from those image patches.
To make full use of the correlations between multicoil images, we propose to introduce the joint sparsity of multiple coils, and we formulate this problem as follows: where w 0,2 = L l=1 |w l | 2 0 is a joint 0 -norm similar to the joint 1 -norm [29], [43].

B. PROBLEMS SOLUTIONS
By introducing the auxiliary variables x = F −1 f ∈ C LQ and B j = W P j (x) ∈ C n×L , and their corresponding Lagrange multipliers (u x ∈ C LQ , and u B j ∈ C n×L , respectively), problem (12) can be transformed into the following subproblems VOLUME 8, 2020 by using the ADMM method: where B ∈ C n×LQ is a matrix concatenating Transform Update: let X = [P 1 (x), . . . , P j (x), . . . , P Q (x)] ∈ R n×LQ ; then, problem (15) can be rewritten as: According to references [11], [13], if X (B − u B ) H has a full singular value decomposition (SVD) of U V H , the solution of problem (20) is given by: Sparse coding: problem (15) can be written as: The solution of problem (22) is given by: where H (x, θ) denotes the hard-thresholding operator [11], which is defined as: where x and θ represent the input matrix and threshold, respectively. When considering the joint sparsity of multiple coils, problem (22) should be modified to the following joint denoising problem: Problem (25) can be solved by: where H J (x, θ) is defined as: Image update: Let the derivative of the objective function of (16) be zero; then, we can obtain: Since W H W = I and Q j=1 P H j P j is a diagonal matrix, the analytical solution is given by: Reference [44] introduces a majorizer for the function J r (R x (f )) at pointf i−1 : where N r (X ) is an operator that generates a 2LN R ×(2LN R − r) matrix whose columns are equal to the right singular vector associated with the smallest q−r value or zero singular value in the generalized SVD of X . Then, problem (17) can be reformulated as follows: Problem (31) can be solved using the CG method. All subproblems can be solved, so we obtain parallel MR imaging reconstruction algorithms based on transform learning and joint transform learning called: PLORAKS-TL and PLORAKS-JTL, which are summarized in Algorithm 1.

IV. SIMULATION EXPERIMENT AND ANALYSIS A. EXPERIMENTAL SETUP
In the following experiments, we compare the PLORAKS-TL and PLORAKS-JTL algorithms with the PLORAKS (with no regularization term), and PLORAKS-JTV (with the JTV regularization term) models. All the compared algorithms are implemented using MATLAB (MathWorks, Natick, MA).
To compare the reconstruction performance of all the algorithms, we used four sets of fully-sampled in vivo human data sets, namely dataset1, dataset2, dataset3, and dataset4, as shown in Fig. 1. Dataset1 and dataset2 are human brain slices acquired by using an eight-channel head coil. Dataset3 is cine data set acquired by using a 28-channel coil,

Algorithm 1 Joint Sparsity and Transform Learning Based PLORAKS Parallel MR Imaging Reconstruction Algorithms (PLORAKS-TL and PLORAKS-JTL)
H J W i+1 P j (x) + u B j , α/µ 2 for PLORAKS-JTL 7: x i+1 is given by using the equation (29) 8: f i+1 is given by solving the problem (31) using the conjugate gradient method 9: 11: until some stopping criteria met 12: Output: f = f i+1 which is then compressed to 12 virtual coils. Dataset4 is a knee data set acquired by an 8-channel HD knee coil. The testing data sets are generated by retrospectively undersampling the fully-sampled data sets using the 2D Poisson-disc undersampling mask with an acceleration factor of R (excluding ACSs), as shown in Fig. 2.
All the following experiments are carried out on a laptop computer with an Intel (R) Core (TM) i5-10210U processor @ 1.6 GHz, 12 GB memory and a Windows 10 operating system (64 bit). All the parameters of the compared algorithms are tuned to achieve the optimal signal-to-noise ratio (SNR) [30] values. In the following experiments, we use the SNR, normalized root mean square error (NRMSE) [28], high-frequency error norm (HFEN) [10], and structural similarity index measure (SSIM) [45] as metrics to quantitatively evaluate the quality of reconstructed images. Higher SNR and SSIM values or lower NRMSE and HFEN values indicate better image reconstruction quality. The SNR is defined as follows: 10 Var MSE (32) where MSE represents the mean square error between the reconstructed imagex and the reference image x, Var represents the variance of x, and x andx are calculated from the fully sampled k-space data and the reconstructed k-space data, respectively by using SoSF(f ). SoSF(f ) is defined as follows: The NRMSE is defined as follows: The HFEN is defined as follows: where LoG (·) is a Laplacian Gaussian filter, which is used to capture the edges of the image, the size of the filter kernel is 15 × 15 pixels, and the standard deviation is 1.5 pixels. The SSIM is defined as follows: where u x and ux are the means of x andx, respectively, σ 2 x and σ 2 x are the variances of x andx, respectively, and σ xx is the covariance ofx and x.    edges in the reconstructed images, and the PLORAKS-TL and PLORAKS-JTL models further alleviate the artifacts and preserve the details missing in the images output by other algorithms. As shown in Fig. 6, the PLORAKS-JTL model solves the problem of oversmoothed edges observed in the image reconstructed by the PLORAKS-TL model.
We quantitatively evaluate the quality of the images reconstructed by the four compared algorithms using 2D  Poisson-disc undersampling with different sampling ratios and 24 × 24 ACSs for dataset1-dataset4. Tables 1-4 tabulate the SNRs, NRMSEs, HFENs, and SSIMs of all the compared algorithms.
As shown in Tables 1-3, PLORAKS-JTV can achieve average SNRs that are approximately 1 dB, 1.4 dB, and 0.3 dB higher than those of PLORAKS for dataset1, dataset2, and dataset3, respectively. The PLORAKS-TL and PLORAKS-JTL models can obtain comparable average SNR values, which are approximately 1.3 dB, 1.8 dB, and 0.9 dB higher than those of the PLORAKS-JTV model, for dataset1, dataset2, and dataset3, respectively. As shown in Table 4, PLORAKS-JTV can achieve an average SNR that is approximately 1.6 dB higher than that of PLORAKS. The PLORAKS-TL model can obtain an SNR that is approximately 2.3 dB higher than that of PLORAKS-JTV. The average SNR of PLORAKS-JTL is approximately 0.9 dB  higher than that of PLORAKS-TL. Overall, introducing an additional regularization term contributes to improving the quality of the image reconstructed by PLORAKS. The PLORAKS-JTL algorithm always achieves the best performance in terms of the SNRs, NRMSEs, HFENs, and SSIMs among the compared algorithms.

C. ALGORITHMIC PERFORMANCE COMPARISON WITH DIFFERENT ACS SIZES
In addition, PLORAKS is a calibration-less parallel image reconstruction method that does not require ACSs [38]. We also compare PLORAKS and PLORAKS-JTL with a coil-by-coil auto-calibrating reconstruction method using the joint TV regularization term (SPIRiT-JTV) [31] from the undersampled k-space data with an acceleration factor of 3 and different ACS sizes. As shown in Table 5, SPIRiT-JTV can achieve an SNR of 21.24 dB when the size of the ACSs   is 24 × 24; when the size of the ACSs is reduced to 12 × 12, the SNR value decreases by approximately 2 dB; and the SNR value dramatically decreases to 10.04 dB when the ACSs is 8 × 8. Therefore, we do not perform reconstruction experiments for the 4 × 4 and 0 × 0 ACS sizes. For the PLORAKS algorithm, when the size of the ACSs gradually reduced from 24×24 to 0×0, the SNR of reconstructed image decreases from 20.85 dB to 17.52 dB, with a reduction of only 3.33 dB. PLORAKS-JTL can effectively improve the quality of images reconstructed by PLORAKS and reduce the impact of the size of the ACSs on the reconstruction performance. Specifically, when the size of the ACSs decreases from 24×24 to 0×0, the SNR of the reconstructed image decreases from 22.23 dB to 20.77 dB, with a reduction of only 1.46 dB.
In general, the SPIRiT-JTV method depends heavily on the size of the ACSs. When the size of the ACSs gradually decreases from 24 × 24 to 0 × 0, the SNR of the image

D. SENSITIVITY TO PARAMETER SETTINGS
In this subsection, we discuss the sensitivity of the proposed algorithm's parameters by adjusting one parameter and fixing the others. The regularization parameter α exists in the form of α/µ 2 in our method, and we set α as 2 × 10 −7 . We have found that β and µ 1 are insensitive to reconstruction accuracy, so we set β as 10 −4 and µ 1 as 10 −4 . We then evaluate the proposed algorithm on dataset1 with a 2D Poisson sampling mask with 4× acceleration and 24 × 24 ACSs using different parameter sets, i.e., µ 2 ∈ [10 −4 , 10 −1 ] and λ ∈ [10 −6 , 10 −3 ], logarithmically spaced. Contour plots show the values of the SNR, NRMSE and HFEN versus λ and µ 2 when reconstructing dataset1 from the Poisson-disc undersampling with 4× acceleration and 24 × 24 ACSs, as shown in Fig. 7. We choose λ and µ 2 that maximize the SNR and minimize the NRMSE and HFEN, respectively. The optimal λ = 10 −5 and the optimal µ 2 = 10 −3 ; and these values are roughly consistent for all three metrics (SNR, NRMSE and HFEN).

E. CONVERGENCE PROPERTY
To reflect the convergence property of the proposed method (PLORAKS-JTL), we explore the development of three metrics (SNR, NRMSE and HFEN) with respect to the number of iterations. Fig. 8 illustrates the values of the SNR, NRMSE and HFEN versus the number of iterations when reconstructing dataset1 from the 2D Poisson-disc undersampling with 4× acceleration and 24 × 24 ACSs. As shown in Fig. 8, all three metrics achieved approximate convergence after 30 iterations.

F. COMPARISON OF RUN TIME
We also compare the run time of PLORAKS, PLORAKS-JTV, PLORAKS-TL, and PLORAKS-JTL. Table 6 tabulates the run time of all the compared algorithms for parallel MR reconstructions from the 2D Poisson-disc undersampling with 6× acceleration and 24 × 24 ACSs for the dataset1. As shown in Table 6, the PLORAKS-JTL method is the slowest among the compared algorithms. The run time of PLORAKS-JTL is K out × (K CG × T CG + T JTL ), where K out is the number of outer iterations, K CG is the number of CG iterations, T CG is the run time of one CG iteration, and T JTL is the run time of one JTL denoising iteration. In the future, we will try our best to solve problem (13) more efficiently, and we will optimize the algorithm by using the graphics processing unit (GPU) to greatly reduce the run time of each JTL denoising iteration.

G. COMPARISON RESULTS WITH A 1D UNDERSAMPLING SCHEME
We also compare the MR image reconstructions for dataset1 from a 1D uniform undersampling and 1D random undersampling with 3× acceleration and 20 ACS lines, as shown in Fig. 9. The visual comparisons are shown in Figs. 10-11. As shown in Figs. 10-11, we can see that PLORAKS-TL and PLORAKS-JTL greatly improve the imaging accuracy and reduce the artifacts exhibited in the images output by other algorithms, and PLORAKS-JTL has slightly better imaging quality than PLORAKS-TL.

V. CONCLUSION
In this paper, we propose two algorithms that combine joint sparsity and adaptive sparsifying transform learning with the PLORAKS model to improve the reconstruction quality of parallel MR imaging, and the reconstruction problems are solved by using the ADMM and conjugate gradient techniques. The experimental results show that the two proposed algorithms can achieve higher performance than the PLO-RAKS algorithm and the PLORAKS-JTV algorithm with the joint TV regularization term in terms of their SNR, NRMSE, HFEN, and SSIM values. In addition, the proposed algorithms are insensitive to the size of the ACSs.