Sparse Constrained Reconstruction for Accelerating Parallel Imaging Based on Variable Splitting Method

Abstract

Parallel imaging is a rapid magnetic resonance imaging technique. For the ill-conditioned problem, noise and aliasing artifacts are amplified during the reconstruction process and are serious especially for high accelerating imaging. In this paper, a sparse constrained reconstruction problem is proposed for parallel imaging, and an effective solution based on the variable splitting method is contrived. First-order and second-order norm optimization problems are first split, and then they are transferred to unconstrained minimization problem by the augmented Lagrangian method. At last, first-order norm and second-order norm optimization problems are alternatively resolved by different methods. With a discrepancy principle as the stopping criterion, analysis of simulated and actual parallel magnetic resonance image reconstruction is presented and discussed. Compared with the routine parallel imaging reconstruction methods, the results show that the noise and aliasing artifacts in the reconstructed image are evidently reduced at large acceleration factors.

1. Introduction

Parallel imaging is a robust method for accelerating the acquisition of magnetic resonance imaging (MRI) data, which exploits spatial sensitivity of an array of receiver coils to reduce the number of the required Fourier encoding steps. However, these reduced amounts of MR data lead to aliased images by the routine reconstruction method. Over the past few years, a number of parallel MRI techniques have been proposed for reconstructing MR image from these undersampled data in either 𝑘-space or image domain [1]. Sensitivity encoding (SENSE) [2] and generalized auto-calibrating partially parallel acquisitions (GRAPPA) [3] are two methods most commonly used on clinical scanners today.

As the amount of data acquired in parallel MRI is less, which depends on the acceleration factor (AF), the quality of the reconstructed image is poorer. Therefore, the AF is usually lower when parallel imaging technique is used to speed up MRI in clinic. Based on SENSE method, the image reconstruction for parallel imaging is performed by solving a linear system that explicitly depends on the sensitivity maps of the receiver coils, and some prior information may be used to stabilize the reconstruction process. Regularization is an attractive means of restoring stability in the reconstruction mechanism, where prior information can be effectively incorporated [4]. The Tikhonov regularization is a commonly used method where a low-resolution prior image is applied in the reconstruction; that is, a quadratic minimization problem and its numerical algorithm are simple, such as the linear conjugate gradient (CG) method. However, the disadvantage of Tikhonov regularization method is that it biases the estimated reconstructed image towards the prior image [5]. More recently, total variation methods have been investigated for MR image reconstruction. The advantage of this type of regularization is that it biases the reconstructed image towards a piecewise smooth image, instead of a globally smooth image, thereby better preserving image edges [6]. With the advent of compressed sensing (CS) theory, sparsity-promoting regularization criteria have gained popularity in MRI, which is known as sparseMRI or CS-MRI [7]. The basic assumption underlying CS-MRI is that many MR images are inherently sparse in some transform domain and then can be reconstructed with high accuracy from significantly undersampled 𝑘-space data. Certainly, the CS framework is apt for pMRI with undersampled MR data [8].

This paper investigates the problem of sparse constrained reconstruction from highly undersampled MR data for parallel MRI. Based on sparseMRI theory, we use the finite difference as the sparsity project domain, and sparsity property of anisotropic total variation (TV) of MR image is used as the prior information for stabilizing the reconstruction process. As a result, a nonlinear optimization problem is constructed for reconstructing the parallel MR image. However, as there exists both first-order norm (known as ℓ1-norm) and second-order norm (ℓ2-norm) minimization problems, the solution to the constructed optimization problem is commonly difficult. In the paper, the ℓ1-norm minimization problem for sparse representation of MR image and ℓ2-norm minimization problem, which are subject to data consistency based on SENSE method, are firstly split by variable splitting method, and then the split constrained minimization problem is converted to an unconstrained minimization problem by the augmented Lagrangian (AL) method. At last, the Lagrangian multiplier method and alternating direction method (ADM) are used to solve the split minimization problem by the different numerical algorithms. In order to evaluate the effectiveness and robustness of the proposed algorithm, the image reconstruction problem from highly undersampled parallel MR data is exploited.

2. Theory

2.1. The MR Signal Model in Parallel Imaging

As a fast imaging method, parallel imaging technique is also known as multiple coils MRI, which uses an array of RF receiver surface coils to acquire multiple sets of undersampled 𝑘-space data simultaneously. Let 𝑟⃗ denote the two-dimension spatial coordinates (𝑥,𝑦) and 𝑠_𝑙(𝑟)⃗ the demodulation information associated with the 𝑙th coil, then the MR signal associated with the 𝑙th coil is the following general forward model: 𝑐_𝑙(𝑟)⃗=𝑠_𝑙(𝑟)⃗𝑓(𝑟)⃗𝑒^{⃗𝑟−𝑖𝑘⋅}𝑑𝑟,⃗𝑙=1,…,𝐿. (1)

Here, 𝑓(𝑟)⃗ denotes the object’s transverse magnetization signal and forms MR image. Let 𝑦-ordinates denote the gradient encoding direction, consideration of MR signal sampled on Cartesian coordinates, and then the discrete MR signal is obtained as follows (2) in accelerating parallel imaging: 𝐼_𝑙𝑥,𝑦=^𝑁𝐴−1_𝑛=0𝑠_𝑙𝑥,𝑦+𝑛𝑀𝑓𝑥,𝑦+𝑛𝑀, (2) where 𝑁_𝐴 is the number of aliased pixels; 𝑅 is known as AF; 𝑀=𝑁_𝑦/𝑅; 𝑥=0,…,𝑁_𝑥; 𝑦=0,…,𝑁_𝑦/𝑅; 𝑁_𝑥, 𝑁_𝑦 is the discrete pixels number, respectively, along the 𝑥-direction and 𝑦-direction when data are full sampled.

2.2. The Model of Image Reconstruction in Parallel Imaging

Considering of noisy samples of MR signal, the discrete model for parallel imaging is given as 𝑦=𝐹𝑆𝑓+𝜀, (3) where 𝑓 is a 𝑁×1 column vector containing the samples of the unknown image to be reconstructed, 𝑦 and ε are 𝑀𝐿×1 column vectors corresponding to the data samples from 𝐿 coils and noise, respectively. 𝑆 is 𝑁𝐿×𝑁 matrix given by 𝑆=[𝑠^𝐻₁,…,𝑠^𝐻_𝐿], 𝑠_𝑙 is 𝑁×𝑁 diagonal matrix corresponding to the sensitivity map of the 𝑙th coil, and 1≤𝑙≤𝐿, (⋅)^𝐻 represents the Hermitian-transpose. 𝐹 is 𝑀𝐿×𝑁𝐿 matrix given by 𝐹=𝐼_𝐿⊗Fu, Fu is 𝑀×𝑁 Fourier encoding matrix, 𝐼_𝐿 is the identity matrix of size 𝐿, and ⊗ denotes the Kronecker product. In order to speed up parallel imaging, 𝑘-space MR data may be undersampled to reduce the total scan time, so 𝑀≤𝑁.

Given an estimate of the sensitivity maps 𝑆, the image reconstruction problem for parallel imaging is to find 𝑓 from data 𝑦. Based on SENSE method, the MR image may be reconstructed using the least-squares estimation as 𝑓=𝑆^𝐻𝐹^𝐻Ψ⁻¹𝐹𝑆⁻¹𝐹^𝐻𝑆^𝐻Ψ⁻¹𝑦, (4) where Ψ denotes the 𝐿×𝐿 receiver noise matrix, it describes the levels and correlation of noise in the receiver channels. Because of the coil configurations and the coil sensitivity error, the measurement matrix FS is commonly non-orthogonal as 𝑆⁻¹𝐹⁻¹𝐹𝑆≠𝐼. The linear system as (3) is ill-posed. Using inversion matrix method as (4), noise derived from measurement data may be amplified if the small eigenvalues exist in matrix 𝐹𝑆, which might result in the instability of reconstruction process. As noise in MRI measurements is Gaussian distribution, a natural approach is to estimate 𝑓 by minimizing a regularized least-squares cost function: 𝑓=argminΦ𝑓=𝐹𝑆𝑓−𝑦₂+𝜇𝑅𝑓, (5) where 𝑅(𝑓) denotes the regularization term, and ‖‖₂ is ℓ2-norm, which represents the data-consistent term. 𝜇 is the so-called regularization parameter to balance the regularization term and data-fidelity term.

An open problem in most regularization image reconstructions is how to best choose the regularization term. If this term is not included, then the image estimate will suffer from noise and aliasing artifacts for undersampled data. The simplest choice is the Tikhonov regularization 𝑅(𝑓)=||𝑓||₂ or 𝑅(𝑓)=||𝑓−𝑓₀||₂, where 𝑓₀ is a prior or reference image. However, if the reference image is zero, then all pixel values in estimation 𝑓 are diminished towards zero, possibly reducing contrast. Another choice is a quadratic roughness penalty function, by which it is convenient for minimization and guarantees that the cost function (5) has a unique minimization, but it has the drawback of smoothing image edges. TV regularization reconstruction may overcome the drawbacks mentioned above, but it is getting harder to be minimized and can lead to the appearance of “blocky” texture in images.

2.3. The Sparse Constrained Image Reconstruction in Parallel Imaging

Based on CS-MRI theory, MR image on finite difference domain is sparse, on which the sparse representation of MR images can be demonstrated by applying a sparsifying transform to a fully sampled image. In this paper, TV transform is used to transform the estimated MR image to finite difference domain, and then the sparse constrained minimization problem for parallel MR image reconstruction is obtainedas min_𝑓∇𝑓₁s.t.𝐹𝑆𝑓−𝑦₂<𝜎, (6) where ‖‖₁ denotes ℓ1-norm, and σ represents any error, such as noisy level derived from the data sampled process. ||∇𝑓||₁ is called TV norm of an image, which is defined as a function of the image gradient. ||∇𝑓||₁ might be computed by discrete isotropic TV or anisotropic TV as ∇𝑓₁=_𝑖∇_𝑥𝑓²_𝑖+∇_𝑦𝑓²_𝑖,∇𝑓₁=_𝑖∇_𝑥𝑓+∇_𝑦𝑓, (7) where ∇_𝑥𝑓 denotes the 𝑥-direction gradient of image 𝑓, and ∇_𝑦𝑓 denotes the 𝑦-direction gradient. There are both ℓ1-norm term and ℓ2-norm term in (6). The ℓ1-norm term in (6) is based on the sparsity of MR image on finite difference domain, while ℓ2-norm term is based on the data fidelity of MR image reconstruction. Using the common numerical algorithms such as CG or Newton method, it is difficult to attain the stabilized solution to problem (6). Lustig et al. researched sparseMRI and posed the detailed computation method for sparse constrained inverse problem [8], where nonlinear CG method was applied to solve the constrained minimization problem as (6) when isotropic TV is applied to compute TV. However, the reconstructed image obtained by this method is dissatisfied according to our analysis.

According to the research product about ℓ1-norm minimization problem [9–11], the variable split method is adopted to split the ℓ1-norm term and ℓ2-norm term in (6) in our research. Specifically, with an auxiliary variable 𝑤, let 𝑤=∇𝑓, then the constrained optimization problem is obtained, which denoted as 𝑃 as follows: 𝑃∶min_𝑤,𝑓‖𝑤‖₁s.t.𝑤=∇𝑓,𝐹𝑆𝑓−𝑦₂<𝜎. (8)

3. The Solution to Sparse Constrained Image Reconstruction Based on Variable Splitting Method for Parallel Imaging

3.1. The Unconstrained Parallel MR Image Reconstruction Problem

In the augmented Lagrangian (AL) framework (also known as the multiplier method [12]), an AL function can be constructed for problem (8) as min_{𝑤,𝑓,𝛾1,𝛾2}Φ𝑓,𝑤,𝛾₁,𝛾₂=‖𝑤‖₁+𝛾^𝐻₁+𝜇𝐹𝑆𝑓−𝑦₁2𝐹𝑆𝑓−𝑦₂+𝛾^𝐻₂+𝜇𝑤−∇𝑓₂2𝑤−∇𝑓₂, (9) where 𝛾₁, 𝛾₂ represent the vector of Lagrange multipliers, and 𝜇₁, 𝜇₂ are the regularization parameters. The solution to (9) may be in the following AL version: 𝑓^𝑘+1,𝑤^𝑘+1⟵argmin_𝑓,𝑤Φ𝑓,𝑤,𝛾^𝑘₁,𝛾^𝑘₂,𝛾₁^𝑘+1⟵𝛾^𝑘₁−𝜆𝜇₁𝐹𝑆𝑓^𝑘+1,𝛾−𝑦_2𝑘+1⟵𝛾^𝑘₂−𝜆𝜇₂𝑤^𝑘+1−∇𝑓^𝑘+1, (10) where 𝜆∈(0,2) guarantees convergence, as long as the subproblem is solved to an increasingly high accuracy at every iteration.

The joint minimization of Φ with respect to 𝑓 and 𝑤 can be computationally challenging in (10). ADM [13] is applied, which alternatively minimizes Φ with respect to one variable at a time while holding others constant. This method decouples the individual updates of 𝑓 and 𝑤 and simplifies the optimization task. Specifically, at the 𝑘th iteration, we perform the following individual minimizations, taking care of using updated variables for subsequent minimizations and the following algorithm: 𝑤^𝑘+1⟵argmin_𝑤Φ𝑓^𝑘,𝑤,𝛾^𝑘₂,𝜇₂,𝑓^𝑘+1⟵argmin_𝑓Φ𝑓,𝑤^𝑘+1,𝛾^𝑘₁,𝜇₁,𝛾^𝑘₂,𝜇₂,𝛾₁^𝑘+1⟵𝛾^𝑘₁−𝜆𝜇₁𝐹𝑆𝑓^𝑘+1,𝛾−𝑦_2𝑘+1⟵𝛾^𝑘₂−𝜆𝜇₂𝑤^𝑘+1−∇𝑓^𝑘+1. (11)

3.2. The Solution to Minimization Problem with Respect to 𝑤

Holding variable 𝑓 constant, we get the minimization problem with respect to 𝑤 in (9) at the 𝑘th iteration as 𝑤^𝑘+1=argmin_𝑤‖𝑤‖₁+𝛾^𝑘₂^𝐻𝑤−∇𝑓^𝑘+𝜇₂2𝑤−∇𝑓^𝑘₂. (12)

Equation (12) is a ℓ1-norm minimization problem, the solution of which can be estimated by shrinkage rule [14] as 𝑤^𝑘+1=Shrink∇𝑓^𝑘−𝛾^𝑘₂𝜇₂,1𝜇₁∗𝜇₂≜max∇𝑓^𝑘−𝛾^𝑘₂𝜇₂−1𝜇₁∗𝜇₂,0⋅sgn∇𝑓^𝑘−𝛾^𝑘₂𝜇₂. (13)

3.3. The Solution to Minimization Problem with Respect to 𝑓

Holding variable 𝑤 constant, we get the minimization problem with respect to𝑓 at the 𝑘th iteration based on (9) as 𝑓^𝑘+1=argmin_𝑓𝛾^𝑘₁^𝐻+𝜇𝐹𝑆𝑓−𝑦¹2𝐹𝑆𝑓−𝑦₂+𝛾^𝑘₂^𝐻𝑤^𝑘+𝜇−∇𝑓²2𝑤^𝑘−∇𝑓₂. (14)

Ignorant of irrelevant constant, (14) can be written as (15) 𝑓^𝑘+1=argmin_𝑓𝜇₁2𝐹𝑆𝑓−𝑦−𝜂₁₂+𝜇₂2𝑤^𝑘−∇𝑓−𝜂₂₂, (15) where 𝜂₁=𝛾₁/𝜇₁, 𝜂₂=𝛾₂/𝜇₂. Equation (15) is a ℓ2-norm minimization problem, which can be solved by nonlinear conjugate gradient (NCG) descent algorithm with backtracking line search algorithm. During the iterative process, the iterative step may be calculated by inexact line search or the Barzilai and Borwein (BB) method [15].

3.4. The Solution to Sparse Constrained Image Reconstruction for Parallel Imaging

Combining the results in Sections 3.1, 3.2, and 3.3, we now present the algorithm for solving the unconstrained optimization problem (9) as follows.(1) Initialize 𝑓⁰, 𝑤⁰ and regularization parameters 𝜇₁, 𝜇₂.(2) Precompute 𝐹^𝐻𝑆^𝐻𝑓, and let 𝜂⁰₁=𝜂⁰₂=0, 𝑘=0.  Repeat the following.(3) Obtain an update 𝑤^𝑘+1 using an appropriate technique as described in Section 3.2.(4) Update 𝑓^𝑘+1 using an appropriate technique as described in Section 3.3 by NCG method, and the inexact line research method is used to calculate the iterative step.(5) Update the Lagrangian multiplier 𝛾₁ by the following function: 𝛾₁^𝑘+1⟵𝛾^𝑘₁−𝜆𝜇₁𝐹𝑆𝑓^𝑘+1−𝑦. (16)(6) Update the Lagrangian multiplier 𝛾₂ by the following function: 𝛾₂^𝑘+1⟵𝛾^𝑘₂−𝜆𝜇₂𝑤^𝑘+1−∇𝑓^𝑘+1. (17)(7) Let 𝑘=𝑘+1.

Until some stop criterion is met.

4. Experiments

4.1. Data Simulation and Acquisition

In all our experiments, we considers 𝑘-space MR data assigned on the Cartesian grid, and the undersampling pattern is uniformity, so Fourier encoding matrix corresponds to an undersampling version of the DFT matrix. Two sets of parallel MR data are applied to test the proposed algorithm. One set of data is acquired from simulating parallel MRI system, and the other is actual parallel MR data.

In order to construct the simulating multichannel MRI receive system, we considered a noise-free 256 × 256 T1-weighted MR image obtained from the BrainWeb database [16] and then downsampled to 128 × 128. The simulating system had four-channel linear phased-array coil wrapped around the whole brain circumferentially. The coil sensitivity maps were calculated using Biot-Savart’s law [17]. Multichannel images were then created by multiplying T1-weighted MR image with the simulated sensitivity profiles. At last, these images were Fourier transformed to generate the multichannel 𝑘-space data whose size was 128×128×4. In the research, complex Gaussian noise would be added to data set for simulating noisy parallel MR data.

The actual fully sampled brain dataset was obtained from PULSAR (a Matlab toolbox for parallel MRI) [18], which was acquired using MR systems with eight-channel head array and multichannel receiver from a healthy male volunteer with fast spoiled gradient-echo sequence, TR/TE = 300/10 ms, matrix size = 256×256, tip angle = 15°, and FOV = 22×22 cm.

To simulate the accelerating imaging procedure for parallel MRI, 𝑘-space dataset along phased encoding direction might be uniformly decimated to produce undersampling datasets. Additionally, central 12 phase-encoding lines were preserved to be used for coil sensitivity estimation [19]. We mainly researched parallel MR image reconstruction problem when the undersampling rate (also known as AF) attained its maximum value that was the number of receive coils at theory. The size of undersampled data in simulated data analysis is 128×32×4, while actual dataset is 256×32×8. Meanwhile, for two sets of data, the full 𝑘-space data was acquired and then sum-of-squares (SoS) reconstruction [20] was used as a reference image.

Figure 1(a) shows the data collected from the simulating parallel MRI, which will be applied to test the proposed algorithm.

(a)

(b)

(c)

Figure 1 

Simulated pMRI. (a) is the T1-weighted brain MR image, (b) is SoS of individual receiver coil sensitivity profiles calculated based on Biot-Savart’s law, and (c) is SoS reconstruction from undersampled data when AF is 4.

4.2. Calibration of Sensitivity Profiles

Calibration of the spatial sensitivity functions of coil arrays is a crucial element in parallel MRI. Accurate coil sensitivity information is required for accurate spatial encoding in parallel MRI reconstructions, and the choice of sensitivity calibration strategy is at least as important as the choice of reconstruction strategy. The most common approach has been to measure coil sensitivities directly using one or more low-resolution images acquired before or after accelerated data acquisition. However, since it is difficult to ensure that the patient and coil array will be in exactly the same positions during both calibration scans and accelerated imaging, this approach can introduce sensitivity miscalibration errors into parallel MRI reconstructions. Coil sensitivity functions vary slowly as a function of spatial position, and low-resolution in vivo images suffice to form sensitivity references. Thus, valid coil sensitivities can be determined from the fully sampled region of the central 𝑘-space, so long as the range of spatial frequencies covered in the central region contains the spatial frequency band of the coil sensitivity functions [19]. In the paper, central 12 phase-encoding lines were Fourier transformed to produce low-resolution 𝑓(𝑟)⃗, which was used to estimate the coil sensitivity profiles. By (18), the magnetization distribution low-resolution 𝑓(𝑟)⃗ may be partly removed, and then the encoding effects of relatively pure coil sensitivities are isolated as 𝑠_l̂(𝑟)⃗≈𝑓(𝑟)⃗𝑠_l(𝑟)⃗^{low-resolution}∑_𝑙𝑓(𝑟)⃗𝑠_l(𝑟)⃗^{low-resolution2}, (18) where the “low-resolution” superscript indicates that use of only the central 𝑘-space positions results in a low-resolution measurement of the full product of 𝑓(𝑟)⃗ and 𝑠_𝑙(𝑟)⃗.

4.3. Data Analysis

In the study of actual parallel MRI, the standard reference image is the SoS reconstructed image of full-sampled individual coil surface images, while the T1-weighted brain image is the reference image in the simulating study. The normalized mean squared error (NMSE) between the reconstructed image and the reference image would be calculated to quantitatively analyze the quality of reconstructed image. NMSE is also known as artifact power (AP), which suggests any error as both increased image artifacts and noise. As shown in (19), a higher value of NMSE (or AP) represents the reduced image quality as ∑NMSE=^⃗𝑟𝑓(𝑟)⃗−𝑓^reference(𝑟)⃗²∑_⃗𝑟𝑓^reference(𝑟)⃗², (19) where 𝑓^reference denotes the sensitivity reference image.

To pose the effectiveness of the proposed reconstruction method as (9), the other parallel MR image reconstruction methods known as generalized encoding matrix (GEM) method [21], which is used in SENSE method, and CG method [22] are comparatively analyzed. Additionally, the stopping criterion for the iterative algorithms used in our study is based on discrepancy principle as norm𝑓_𝑗−𝑓_𝑗−1norm𝑓_𝑗<tol, (20) where description “norm” indicates Frobenius norm of matrix, “tol” denotes the convergence tolerance, and subscription “𝑗” indicates the iteration number during the iterative process.

4.4. Results

We compared the proposed method to GEM method, named the basic SENSE, and to CG method, by which MR images may be reconstructed from undersampled MR data. Figure 2 shows the reconstructed images from noise-free simulating data when employing 4-fold accelerating MRI, when AF is 4. Figure 3 shows the reconstructed images from noisy 4-fold undersampled MR data, where the added noise is zero-mean Gaussian noise. In Figures 2 and 3, the three images on the top are the reconstructed images by the methods as mentioned above, and the images on the bottom are the difference profiles between the top three reconstructed images and reference image.

(a) Reconstruction Image by GEM

(b) Reconstruction Image based on CG

(c) Sparse Constrained Reconstruction Image

(d) Difference profile between reconstructed reference image

(e) Difference profile

(f) Difference profile

Figure 2 

Reconstructed images from noise-free undersampling MR data for the simulated parallel MRI, where the undersampling rate is 4, which is the maximum of AF at theory. The top three images are, respectively, reconstructed by GEM method, CG method, and the proposed method as (9), which is solved by the algorithm referred to in Section 3.4. The bottom three images are the difference profiles between the above three reconstructed images and the reference image.

(a) Reconstruction Image by GEM

(b) Reconstruction Image by CG

(c) Sparse Constrained Reconstruction Image

(d) Difference profile

(e) Difference profile

(f) Difference profile

Figure 3 

Reconstructed images from noisy undersampled MR data for the simulated parallel MRI, where the undersampling rate (named AF) is 4.

During the iterative process of reconstructing the images as seen in Figures 2 and 3, we let the convergence tolerance 0.0001. When the proposed algorithm (described in Section 3.4) is used to reconstruct the image in Figure 2, the penalty parameters 𝜇₁, 𝜇₂ are 2⁻²⁷, 2⁻³³, respectively. The number of iterations is 20 for resolving the minimization problem with respect to 𝑓 by the NCG method referred to in Section 3.3, and the total number of iterations for the proposed algorithm is 411, namely, the valuable 𝑘 referred to in Section 3.4. At the same time, the proposed algorithm is used to reconstruct the images in Figure 3, and the penalty parameters 𝜇₁, 𝜇₂ are 2⁻²⁶, 2⁻²², respectively, and the number of iteration is 20 for resolving the minimization problem with respect to 𝑓 by NCG method, and the total number of iterations for the algorithm proposed is 198. Tables 1 and 2 are the mean of NMSE, respectively, in Figures 2 and 3.

Table 1 

Mean of NMSE of reconstruction image (in Figure 2).

Table 2 

Mean of NMSE of reconstruction image (in Figure 3).

As shown in Figure 2, the proposed algorithm may obviously restrain the aliasing artifacts resulted from undersampled MR data. As seen in Figure 3, noise derived from measurement procedure is suppressed and result in the convergence of iterative process by our proposed algorithm.

In analysis of actual parallel MR image reconstruction, AF is set as 8, which is the maximum at theory for parallel MRI system with 8-channel coil array. Figure 4 shows the reconstructed image and the difference profiles between reconstructed images and the reference image, and the mean of NMSE is shown in Table 3.

Table 3 

Mean of NMSE of reconstructions image from actual pMRI data (AF = 8, in Figure 4).

(a) Reference Image

(b) Reconstruction Image by GEM

(c) Reconstruction Image by CG method

(d) Sparse Constrained Reconstruction Image

(e) Aliasing Image

(f) Difference profile

(g) Difference profile

(h) Difference profile

Figure 4 

Reconstructed images from noisy undersampling MR data for in vivo parallel MRI, where the undersampling rate is 8, which is the maximum at theory. (a) is the reference image (SoS reconstruction of full-sampled coil images). (b), (c), and (d) are images reconstructed, respectively, by GEM method, CG method, and the proposed method proposed as (9). (e) is the reconstructed image by SoS method. (f), (g), and (h) are the difference profiles, respectively, between (b), (c), (d), and (a).

In Figure 4, the left image on the top penal is the reference image, and the left image on the bottom is the reconstructed image by SoS method from 8-fold undersampled parallel MR data. During iterative process of reconstructing the images in Figure 4, the convergence tolerance is 0.001. When the proposed algorithm is used to reconstruct the image in Figure 4, the penalty parameters 𝜇₁, 𝜇₂ are 2⁻²⁰, 2⁻¹¹, respectively. The iterative number is 1 for resolving the minimization problem with respect to 𝑓 by NCG method, and the total iterative number of algorithm proposed is 37.

As shown in Figures 3 and 4, the proposed algorithm can effectively reduce the aliasing artifacts and noise information.

5. Discussion

The quality of the reconstructed image for parallel imaging may be affected by the undersampling rate, noise level, the exact estimation of receiver coils sensitivity, 𝑘-space data trajectory, and reconstruction method. We only researched the image reconstruction from undersampled 𝑘-space data acquired along the Cartesian trajectory. However, the proposed algorithm is also fit to reconstruct MR image for arbitrary trajectories.

There may be two approaches to improve the quality of the reconstructed image based on SENSE method: one approach is to improve the accuracy estimation of sensitivity, another approach is to propose the effective reconstruction algorithm. However, when AF is large, only by improving the reliability of the coil sensitivity maps, the suppression of aliasing artifact and noise in the reconstructed image is indistinctive. In order to remarkably improve the quality of the reconstructed image, the constrained reconstruction is generally shown to be an effective method. The goal of our proposed method is to improve the quality of the reconstructed image, when MR data are highly undersampled in parallel imaging. As the coil sensitivity profiles are already obtained, the reconstruction problem in parallel imaging is formulated as solving a set of linear equations based on (3). These equations can be very ill conditioned depending on the coil configurations and sampling trajectories, which further deteriorate the quality of the reconstructed image, especially when high undersampling rate is used for accelerating parallel imaging procedure. Therefore, the ill-conditioned problem in parallel imaging has been partially addressed by optimizing coil geometry, optimizing sampling trajectory, or introducing regularizations. In this paper, the ill-conditioned problem is partly reformed by adding the constrained condition to stabilize the reconstruction process. As anywhere noise exits, the reconstruction method for parallel imaging should be proposed on full consideration of both noise and aliasing artifacts suppressed problem. Based on sparseMRI theory, as MR image is sparse on finite difference domain, and a sparse constrained reconstruction problem for parallel imaging is constructed in the paper. Through the methods proposed and the experiment presented, the following aspects could be taken into account.

(1) Compare the proposed reconstruction method with the other methods for parallel imaging. According to our analysis results, the GEM algorithm for usual SENSE method is usually suitable for noise-free data and the accelerating imaging procedure at low AF. CG method can effectively restrain aliasing artifacts in the reconstructed image and be fast convergence on condition of noise-free undersampled data. However, when MR image may be reconstructed from the noisy and undersampled MR data, NMSE value between the reconstructed image and reference image first rapidly decreases and then gradually increases along with CG iterative reconstruction process. As shown in Figure 3, the least mean of NMSE between the reconstructed image by CG method and the reference image was 0.0762 after 21 numbers of CG iteration, but reached the convergence condition after 975 numbers of iteration. Similarly in Figure 4, the least mean of NMSE of the reconstructed image by CG method was 0.0623, but reached the convergence condition after 1020 number of iteration, and the least mean of NMSE was shown in Table 3.

(2) Reconstruct MR image under the nonnegative constraint condition. In our proposed algorithm, the solution to subproblem with respect to 𝑓 as (15) is a quadratic minimization problem, which may be the following subproblem on condition of nonnegative constraint: min_𝑓𝜇₁2𝐹𝑆𝑓−𝑦−𝜂₁₂+𝜇₂2𝑤^𝑘−∇𝑓−𝜂₂₂s.t.real𝑓>0, (21) where real () indicators the operator of getting the real of complex 𝑓.

Figure 5 shows the iteration process during reconstructing MR image 𝑓 in analysis of actual parallel MRI, where (15) and (21) are, respectively, used in the proposed algorithm described in Section 3.4. As shown in Figure 5, when reconstructing image under the nonnegative constraint condition, the iteration process described in Section 3.4 might rapidly converge, and the mean of NMSE was reduced between the reconstructed image and the reference image.

Figure 5 

Iteration convergence process of reconstructing MR image for actual parallel MRI. The reconstruction method is the proposed algorithm described in Section 3.4, and (15) and (21) are, respectively, used for the subproblem with respect to 𝑓, where NCG descent algorithm with backtracking line search algorithm is used to solve (15) and (21).

(3) Reconstruct MR image at the different AFs for parallel imaging. Parallel MRI utilizes a radio frequency (RF) coil array to simultaneously acquire data from multiple receivers, and acceleration is achieved by a reduced phase encoding 𝑘-space trajectory. The nature of the subsampled 𝑘-space data requires the use of a reconstruction algorithm to restore aliased images into full field-of-view (FOV) images. Table 4 is the mean NMSE of reconstructed image, respectively, by GEM method, CG method and the proposed algorithm at the different AFs.

Table 4 

The mean of NMSE of reconstructed image from actual pMRI data (AF = 4).

As shown in Table 4, the proposed algorithm obviously improves the quality of reconstructed images. Compared to Table 3, AF is the main reason for affecting the quality of reconstructed image, so aliasing artifacts suppression may be chiefly considered for reconstruction method.

(4) The solution to the quadratic minimization problem is as follows. The solution to subproblem with respect to 𝑓 as (15) is a quadratic minimization problem, which can be solved by the steepest descent method, the Newton method, or CG method. However, when resolving the subproblem with respect to 𝑓, theNewton method is inferior to the nonlinear CG method by the data analysis results.

(5) Penalty parameters and iteration convergence problem for regularization reconstruction are as follows. Convergence problem should be firstly considered for any iteration algorithm. The proposed algorithm for reconstructing MR image is an iterative solution course, and penalty parameters 𝜇₁, 𝜇₂ greatly affect the iterative convergence process. As long as penalty parameters chosen are suitable, the proposed algorithm can rapidly converge. Thereby, the penalty parameters 𝜇₁, 𝜇₂ are the key factors of optimization problem of sparse constrained reconstruction as shown in (9), which is a ℓ1-norm and ℓ2-norm minimization problem, named ℓ1-ℓ2 optimization problem. The common L-curve method is unable to determine the penalty parameters 𝜇₁, 𝜇₂. According to the quantitative index as NMSE between the reconstructed image and the reference image, we explored the better suitable penalty parameters within a determinate range. Specifically, the ratio of μ₁ to μ₂ was fixed, and μ₁, μ₂ increased by power of 2 in the range of 2⁻³⁰ to 2¹⁰. The penalty parameters μ₁, μ₂ would be selected in our research, when they were used to (9) and then the least mean of NMSE between the reconstructed image and the reference image can be obtained.

6. Conclusion

In order to improve the quality of the reconstructed image for parallel imaging, a sparse constrained image reconstruction algorithm based on variable splitting method is proposed. Through the analysis of reconstructing full-FOV image from simulating and actual parallel MR data, the proposed algorithm can effectively suppress the aliasing artifacts in the reconstructed image resulted from the undersampled MR data, and also noise is obviously suppressed and the edge of image is preserved.