Joint sparse reconstruction of multi-contrast MRI images with graph based redundant wavelet transform

Background Multi-contrast images in magnetic resonance imaging (MRI) provide abundant contrast information reflecting the characteristics of the internal tissues of human bodies, and thus have been widely utilized in clinical diagnosis. However, long acquisition time limits the application of multi-contrast MRI. One efficient way to accelerate data acquisition is to under-sample the k-space data and then reconstruct images with sparsity constraint. However, images are compromised at high acceleration factor if images are reconstructed individually. We aim to improve the images with a jointly sparse reconstruction and Graph-based redundant wavelet transform (GBRWT). Methods First, a sparsifying transform, GBRWT, is trained to reflect the similarity of tissue structures in multi-contrast images. Second, joint multi-contrast image reconstruction is formulated as a ℓ2, 1 norm optimization problem under GBRWT representations. Third, the optimization problem is numerically solved using a derived alternating direction method. Results Experimental results in synthetic and in vivo MRI data demonstrate that the proposed joint reconstruction method can achieve lower reconstruction errors and better preserve image structures than the compared joint reconstruction methods. Besides, the proposed method outperforms single image reconstruction with joint sparsity constraint of multi-contrast images. Conclusions The proposed method explores the joint sparsity of multi-contrast MRI images under graph-based redundant wavelet transform and realizes joint sparse reconstruction of multi-contrast images. Experiment demonstrate that the proposed method outperforms the compared joint reconstruction methods as well as individual reconstructions. With this high quality image reconstruction method, it is possible to achieve the high acceleration factors by exploring the complementary information provided by multi-contrast MRI.


Background
Multi-contrast images in magnetic resonance imaging (MRI) provide abundant contrast information reflecting the characteristics of the internal tissues of human bodies, and thus have been utilized in clinical diagnosis. However, long acquisition time limits the application of multicontrast MR imaging.
Under-sampling the k-space data and reconstructing images with sparsity constraint is one efficient way to accelerate MRI sampling [1][2][3][4][5]. However, the data acquisition factor is limited since images are compromised when images are reconstructed individually. The previous work [6] suggested to use another fully-sampled contrast image to train an adaptive sparse representation with Graph-based redundant wavelet transform (GBRWT) and then greatly improve the reconstructed images [7]. This approach, however, cannot reduce the overall accel-eration factor in data acquisition because of the full sampling in another contrast images [6]. Thus, to further accelerate multi-contrast MRI, under-sampling all multicontrast images, e.g.T1 weighted (T1W), T2 weighted (T2W) and proton density weighted (PDW) MRI images, and maintain high quality image reconstruction are expected.
The MRI image structures under different contrast settings are the same due to the multiple acquisitions of the same anatomical cross section [6,[8][9][10][11][12]. Thus, non-zero coefficients may occur at the same spatial locations in the sparsifying transform domains, e.g. finite difference, wavelet transform [2] and patch-based sparse transformations [13][14][15][16]. Therefore, it is possible to improve the image reconstruction if this extra information is incorporated into sparse image reconstruction [17].
Sparse representation capability plays a key role in sparse MRI reconstruction. The GBRWT [6,7] transform was verified to have good sparsification capability for MRI images. The main step of GBRWT transform is to construct a graph to find new Fig. 1 Flowchart of the proposed method. The x ref denotes the reference image used to train the graph wavelet transform, y 1 , y 2 , ⋯, y T denote the under-sampled k-space data of multi-contrast images,x 1 ;x 2 ; ⋯;x T denote the reconstructed images Fig. 2 The permutation and wavelet filtering on re-ordered image pixels, a and b are the forward and inverse transforms permutations adaptive to target image structures, and then to obtain the sparser transformation with wavelet filters acting on the permutated smooth signals. However, if high acceleration factor is set, very limited information will be provided for single image thus the reconstruction will be compromised. Thus, the combining merits of joint reconstruction and GBRWT are expected.
In this study, we propose to reconstruct the multicontrast MRI with adaptive GBRWT sparse representations and joint sparsity among multi-contrast images. An alternating direction method with continuity (ADMC) [18] algorithm is introduced to solve the joint ℓ 2, 1 -norm minimization problem. The proposed approach will be compared with the joint sparse reconstruction method based on shift invariant discrete wavelet transform (SIDWT) [17] and Bayesian compressed sensing (BCS) [19].

Methods
The under-sampled k-space data of multi-contrast MRI images are expressed as where x = [x 1 ; ⋯; x T ] denotes the column stacked multicontrast images, T the number of contrasts. y = [y 1 ; ⋯; y T ] the column stacked under-sampled k-space data, ε noises in the sampled data. The UF is the under-sampling operator in the Fourier space, which can be expressed as Each U i F i , (i = 1, ⋯, T) acts on one of the multicontrast images. We adopt different sampling patterns, i. e. U 1 ≠ ⋯ ≠ U i ≠ ⋯ ≠ U T , and the same Fourier transform bases, i.e.
The flowchart of the proposed joint sparse reconstruction is shown in Fig. 1. Reconstructed image based on SIDWT [17] is adopted as reference image to train the GBRWT from the under-sampled k-space data, because SIDWT can mitigate the blocky artifacts  Fig. 4 The grouping operator G and the grouped GBRWT coefficients introduced by orthogonal wavelet transform and better preserve the structures in the target images [15,16]. With GBRWT as the sparse representation, multi-contrast images can be simultaneously reconstructed by implementing joint sparsity constraints on these transformation coefficients.

Graph-based redundant wavelet transform
Given a reference image, the GBRWT is achieved by carrying out redundant wavelet transform on permuted signals of new orders [7]. The new orders are found in weighted graph constructed from the reference image, in which image patches collected by a sliding window serve as the vertex and the patch similarities computed using w m, n = w(b m , b n ) = ‖b m − b n ‖ 2 (where b m and b n denote the m th and the n th patches) serve as the weight. The new orders are obtained by finding the shortest possible path on the patch-based graph [7,20]. Then, redundant wavelet transform is performed on permuted pixels to achieve sparse representation. The process of permutation and wavelet filtering in GBRWT is shown in Fig. 2. In the l th level decomposition, the input signal a l will be first reordered by permutation matrix P l , whose inverse process is P H l and satisfied P H l P l ¼ I . Then, nondecimated wavelet transformation Φ l , whose inverse process is Φ H l and satisfied Φ H l Φ l ¼ I, are performed on the re-ordered pixels. The output a l + 1 and d l + 1 of l th level non-decimated decomposition will be of the same size with the input signalã l . Let Φ l P l be the l th level decomposition of GBRWT, and P H l Φ H l be corresponding composition process, and then Φ l P l satisfies the following property where c denotes the redundancy of GBRWT transform. It has been verified that GBRWT provides sparser representations than traditional wavelet transform, thus can improve the MRI image reconstruction [7].

Joint sparsity of multi-contrast image coefficients
Multi-contrast MRI images are obtained by different parameter settings, but share the same anatomical cross section [6][7][8]. The image structures corresponding to tissue locations remain unchanged with contrast varied, which leads to spatial position-related coefficients. Joint sparsity means that, under appropriate sparsifying transforms, the positions of non-zero coefficients correspond directly to same spatial locations in multiple images. Figure 3 shows that the nonzero transform coefficients of two contrast images occur at the same positions in the Haar wavelets transform and GBRWT domains. Thus, the joint sparsity of GBRWT provides extra information on images and may further improve the reconstruction of multi-contrast images.

Problem formulation
The joint sparsity promoting problem in multi-contrast MRI image reconstruction with GBRWT is solved using the mixed ℓ 2, 1 norm minimization [9,21,22]: where, Ψ ¼ 5 and Ψ g denotes the GBRWT representation, in which l th level decomposition be expressed as Φ l P l . Let α = Ψx be the corresponding coefficients, then for an image set which includes T kinds of contrasts, the column stacked coefficients can be expressed as: α = [α 1 ; ⋯; α T ]. The role of grouping operator G is to reshape the column stacked coefficients of multi-contrast MRI images into a matrix as shown in Fig. 4. Then, one column of Gα stands for coefficients of one image, and one row forms a group.
The ℓ 2, 1 norm is defined as where, G is the group operator, N is the number of transform coefficient and T the number of contrast.

Numerical algorithm
The alternating direction method with continuation [18] is incorporated in the ℓ 2, 1 norm optimization. Let α = Ψx, the objective in Eq. (4) can be rewritten as Furthermore, the objective function in (6) can be overrelaxed to be unconstraint as The λ is a parameter to balance the sparsity and data fidelity. The β is fixed in the inner loops and changes continuously to achieve optimal reconstruction in the outer loops. When β → ∞, the solution of Eq. (7) approaches to that of Eq. (6). When β is fixed, x and α will be computed alternatively by the following two steps: 1) With x fixed, α will be computed by solve the objective: Algorithm 1: Joint multi-contrast MRI reconstruction based on GBRWT Parameters: λ Input: k-space data y = [y 1 ; ⋯; y T ]; g levels of permutation orders P j , {j = 1, ⋯, g}; regularization parameter λ; tolerance of inner loop η = 10 −4 .

Main:
While β ≤ 2 12 (1) Given x, solving α by computing Eq. (10) for each group of coefficients α i , {i = 1, ⋯, N}; (2) Applying α into Eq. (12) to obtain the solution x; (3) If ‖Δx‖ = ‖x previous − x‖ > η, then x previous = x, go to step (1); Otherwise: go to step (4); (4)x¼x, β = 2β, go to step (1) Then, α can be obtained via each group computing: 2) Fix α, x can be computed by solving The minimization with respect to x can be solved by finding the extreme of least square problem in Eq. (11). Finally, we get where, c is the redundancy caused by GBRWT transform. The numerical algorithm pseudo-code is listed in Algorithm 1.

Results
The image reconstruction was performed on a server with E5-2637 v3 (3.5G Hz) *2 CPU, 8 GB memory. The proposed method, Joint sparse reconstruction based on GBRWT (JGBRWT), is compared with the Joint sparse reconstruction method based on SIDWT (JSIDWT) [15][16][17], that replacing Ψ g with SIDWT in Eq. (4) and Joint reconstruction with Bayesian Compressed Sensing (JBCS) [19], which is a state-of-the-art joint multicontrast image reconstruction that jointly explores the gradient coefficients of multiple images. The comparison with GBRWT-based single image reconstruction [7] is also included to demonstrate the advantage of the joint reconstruction. The parameter values for JBCS are taken as the same in the cods shared by the authors (http:// martinos.org/~berkin/software.html). For the proposed method and JSIDWT, λ is set as 10 4 . The relative ℓ 2 norm error (RLNE) defined as eðx _ Þ¼ k x _ −xk 2 =kxk 2 (in whichx is ground truth and x _ is the reconstructed image) and mean structure similarity index Fig. 9 Reconstruction errors of knee images ( × 5 ) with 32% sampled data (Cartesian). Rows 1-3 correspond to multi-contrast images shown in Fig. 6(a-c) respectively. Column 1-4 denote the results obtained from GBRWT, JSIDWT, JBCS and JGBRWT reconstruction respectively measure (MSSIM) [23] served as the criteria for assessing the quality of reconstructed image quality. Smaller RLNE means lower reconstructed error and higher MSSIM indicates better structure preservation capability. The Brainweb images (http://brainweb.bic.mni. mcgill.ca/) [24,25] (Fig. 5) as well as the in vivo multi-contrast images were used to validate the efficiency of the proposed method. The multi-contrast knee images (Fig. 6) were acquired from GE 3 Tesla scanner (Discovery MR750W, USA) with parameters (T1W: FSE, TR/TE = 499 ms/9.63 ms; T2W: TR/TE = 2435 ms/49.98 ms, Proton density weighted image: TR/TE = 2253 ms/31.81 ms; FOV = 180 × 180 mm 2 , slice thickness = 4 mm). The multi-contrast brain images (Fig. 7) were acquired from SIMENS 3 Tesla scanner (MAGNETOM Trio Tim, Germany) with parameters (T2W: TSE, TR/TE = 3000 ms/66 ms,; T1W, FLAIR: TR/TE = 3900 ms/9.3 ms,; FOV = 200 × 200 mm 2 , slice thickness = 5 mm) for Fig. 7(A) and TSE: TR = 4000 ms, FOV = 192 × 192 mm 2 , slice thickness = 3 mm, ΔTE =8 ms for the multi-echo data in Fig. 7(B).  One typical brain image reconstruction with 27% sampled data are shown in Fig. 10. In the zoomedin area of the 2nd row, the sulcus of the T2W image appears in the middle of the fully sampled and JGBRWT reconstructed images, but nearly disappears for JBCS reconstruction. In the marked region of the 3rd, the proposed method leads to more consistent reconstruction with the fully sampled image than other methods. These improvements are also confirmed by the error images in the last two rows.

2D under-sampling
The 2Dunder-sampling patterns (Fig. 11) was explored to demonstrate the potential applications of the proposed method in 3D imaging, in which 2D phase encoding plane can be under-sampled.
Brainweb reconstructed errors shown in Fig. 12 demonstrate that on the simulated database, the lowest reconstruction errors were obtained with the proposed method. The corresponding RLNE/ MSSIM are shown in Table 4. Figure 13 implies that the proposed method led to the lowest brightness in the error images and thus maintained fidelity best. The criteria listed in Tables 4 and 5 indicate that the proposed method achieved the highest MSSIM and the lowest RLNE on the tested dataset.

Different sampling rates
The curves in Fig. 14 show that the RLNEs decreased with sampling rate increased. The RLNE line of the proposed JGBRWT method (dark green line) is lower than that of GBRWT (or contrast-bycontrast reconstruction, black line) with the same GBRWT representation, indicating benefits are achieved by utilizing joint sparsity among multicontrast images. The JGBRWT also outperforms other joint reconstruction method, including JSIDWT (red line) and JBCS (blue line), in terms of lower RLNEs at all sampling rates.

The same sampling patterns
The proposed method is compatible to same or different sampling patterns. Reconstruction criteria in Table 6 show that the proposed method outperforms the compared ones under the same sampling patterns. Besides, at the same sampling rate, using different sampling patterns lead to better evaluation criteria than using same sampling patterns (Table 6 vs. Table 1).

Limitations on choosing image to train the graph
Choose an arbitrary pre-reconstructed image as reference will lead to reconstruction errors (RLNEs) slightly change as shown in Fig. 15. But the RLNEs are still much lower than single image reconstruction. A possible way in the future work is to train a GBRWT jointly from     . 11 The 2D under-sampling mask. a is pseudo-radial sampling with sampling rate 11%, b is random sampling with sampling rate 15% all the under-sampled multi-contrast images to make full use of the common/complementary information of multi-contrast images.

Limitations on un-registered images
Un-registered multi-contrast images will go against the joint sparsity assumption, and thus affect joint reconstruction performance. Reconstructed images of aligned and misaligned multi-contrast images (we simulate misalignment by rotating Fig. 16(a) with 10 degrees) shown in Fig. 16 demonstrate that misalignment will make the detail reconstruction Table 4 RLNE/MSSIM for the reconstruction using 11% pseudoradial k-space data of Fig Reconstructed error ( × 5 ) using 2D Random under-sampling, with 15% sampled data in Fig. 7(a1-a2), and 26% sampled data in Fig. 7(b1-b4). Rows 1-6 are reconstruction errors correspond to multi-contrast images shown in Figs. 7(a1-a2, b1-b4) respectively. Columns 1-4 denote the results obtained from GBRWT, JSIDWT, JBCS and JGBRWT respectively deteriorated. RLNE obviously increased in sparse reconstruction of misaligned multi-contrast images. Improved image reconstruction is expected by incorporating the registration into image reconstruction process as it was done in [6], which would be interesting as a future work.

Computation complexity
The main step of numerical algorithm to solve the proposed joint reconstruction problem include a soft thresholding to solve α and a one-step computation to solve x, which is with the same computation complexity as single contrast image reconstruction, but with more data to compute, and thus no obvious additional computational burden. Program at our platform (E5-2637 v3 (3.5G Hz) *2 CPU, 8 GB memory) shows that, the SIDWT-based single image reconstruction need 20 s, and SIDWT-based joint reconstruction need 100 with 4 different contrast images at low sampling rate. The GBRWT-based single image reconstruction need 200 s and GBRWT-based joint reconstruction need 103 s with 4 different contrast images at low sampling rate.

Experiment with noise
Multi-contrast images in Fig. 7(A) in the manuscript are used in noise experiment. Noisy data are simulated by adding Gaussian white noise with variance σ 2 = 0.02 on real and imaginary part of kspace data. Figure 17 demonstrate that the proposed method outperforms the compared ones in preserving image structures as well as removing noise. According to Table 7 the proposed method achieves lowest RLNEs, highest MSSIMs and highest SNRs. The signal to noise rate (SNR) is defined as SNR = 10log 10 (μ/ σ), where u is the mean of image density and δ is the standard deviation of the noise extracted from the image background.

Parameters
Two noise level are considered (Gaussian white noise with variance σ 2 = 0.02 and σ 2 = 0.03) in testing λ. The optimal λ for σ 2 = 0.02 and σ 2 = 0.03 are 600 and 400 respectively on the tested data according the curve shown in Fig. 18.
The parameters of GBRWT include patch-size and decomposition levels, which have been discussed in [7]. The suggested patch-size in GBRWT are from 4 × 4 to 7 × 7, and suggested decomposition level is 3-5 level. We use the patch-size 7 × 7 and do 5 level decomposition in this experiment. Table 5 RLNE/MSSIM for the reconstruction using 15% randomly sampled k-space data of Fig. 7(A) and Fig. 7