MR Image Reconstruction Based on Iterative Split Bregman Algorithm and Nonlocal Total Variation

Abstract

This paper introduces an efficient algorithm for magnetic resonance (MR) image reconstruction. The proposed method minimizes a linear combination of nonlocal total variation and least-square data-fitting term to reconstruct the MR images from undersampled 𝑘-space data. The nonlocal total variation is taken as the 𝐿₁-regularization functional and solved using Split Bregman iteration. The proposed algorithm is compared with previous methods in terms of the reconstruction accuracy and computational complexity. The comparison results demonstrate the superiority of the proposed algorithm for compressed MR image reconstruction.

1. Introduction

Magnetic resonance (MR) imaging has been utilized in diagnosis because of its glorious depiction of soft tissue changes and noninvasive nature. As explained in [1, 2], it is possible to accurately reconstruct the MR images from undersampled 𝑘-space data using compressed sensing and considerably scale back the scanning period. Suppose 𝑢∈𝑅^𝑛 is a sparse signal, and let 𝐾 be a measurement matrix such that 𝐾𝑢=𝑣, where 𝑣 is the observed data. Then, recovering 𝑢 from 𝑣 is equivalent to solve min_𝑢𝐽(𝑢)suchthat𝐾𝑢=𝑣, (1) where 𝐽(𝑢) is a regularizing function; usually it may be bounded variation or Besov norm. In case of compressed sensing MRI, 𝐾 is partial Fourier matrix (𝐾=𝑃ℱ), where 𝑃∈𝑅^𝑀×𝑁 is an identity matrix (𝑀≪𝑁), ℱ is a discrete Fourier matrix, and 𝑣 is the observed 𝑘-space data contaminated with Gaussian noise of variance 𝜎², the relaxation form for (1) should be given by min_𝑢𝐽(𝑢)suchthat‖𝐾𝑢−𝑣‖²₂<𝜎². (2) In order to make (1) simple to solve, first convert it into an unconstrained optimization problem. One common method for doing this is to use a penalty function/continuation method, which approximates (1) by a problem of the form min_𝑢𝐽(𝑢)+𝜇2‖𝐾𝑢−𝑣‖²₂, (3) where 𝜇 is a balancing parameter between the fidelity term and sparsity term [3, 4]. This is the unconstrained problem we need to solve.

There are a lot of iterative methods existing to reconstruct MR images from undersampled data such as conjugate gradient algorithm (CG) [4], operator splitting algorithm (TVCMRI) [5], variable splitting method (RecPF) [6], composite splitting algorithm (CSA) and its accelerated version called a fast composite splitting algorithm (FCSA) [7], and split Bregman algorithm combined with total variation (SB-TV) [8]. The Split Bregman method provides better solution to a wide class of 𝐿₁-regularized problems. In SB-TV, the Split Bregman technique is applied to the Rudin-Osher-Fatemi functional to a compressed sensing problem that arises in a magnetic resonance imaging. A detailed explanation about the Split Bregman technique is given in Section 2. Reconstruction using CG is very slow for MR images. TVCMRI and RecPF can replace iterative linear solvers with Fourier domain computations, which can provide a reduction in computation time. FCSA performs much better than TVCMRI and RecPF and CSA. FCSA is based on wavelet sparsity and total variation (TV). However, despite the huge popularity of TV minimization, it has also some unwanted effects. In the presence of noise, it tends to piecewise constant solutions, the so called staircasing effect which was analyzed in detail in [9, 10]. Another deficit of TV regularization is the systematic loss of contrast in the reconstructions, even if the given data is noise free. This effect is the well-known systematic error of the total variation minimization and was studied extensively in [11, 12]. The major problem associated with TV-based compressed sensing method is that it tries to uniformly penalize the image gradient irrespective of the underlying image structures, and thus low contrast regions are sometimes over smoothed [13]. To resolve this issue, we propose a new algorithm which jointly minimizes a linear combination of nonlocal total variation and least-square data-fitting term to reconstruct the MR images from under-sampled data. In medical image reconstruction, fine structures, details, and texture should be preserved. The nonlocal total variation makes the recovered image quality sharper, and they preserve the fine structures, details, and textures. NLTV is much better than TV for improving the signal-to-noise ratio in practical application [14–16]. Authors of [17, 18] showed that among the existing nonlocal regularization techniques, NLTV is performing well in reconstruction. So in the proposed method, the nonlocal total variation is taken as the 𝐿₁-regularization functional and solved using Split Bregman algorithm. Numerous experiments have been done on real MR images to show the efficiency and advantages of the proposed work in terms of computational complexity and reconstruction accuracy.

2. Split Bregman Algorithm and Regularization Method

The 𝐿₁-regularized problems have many important applications in engineering and imaging science. The general form of such problems is min_𝑢|Φ(𝑢)|+𝐿(𝑢), (4) where |⋅| represents the 𝐿₁-norm and both Φ(𝑢) and 𝐿(𝑢) are convex functions. Many important problems in image processing can be interpreted as 𝐿₁-regularized optimization problems. Equation (3) is an example of such a problem min_𝑢𝐽(𝑢)+𝜇2‖𝐾𝑢−𝑣‖²₂. (5) Unfortunately, for several issues, selecting a large value for 𝜇 makes (5) extremely tough to solve numerically [19, 20]. Also, for several applications, 𝜇 should be increased in small steps, creating the strategy less efficient. Bregman iteration can be used to reduce (5) into small sequence of unconstrained problems for further processing.

2.1. Bregman Iteration

Bregman iteration is a method for finding extrema of convex functionals [21]. Bregman iteration was already applied to solve the basis pursuit problem in [22] and medical imaging problem in [23]. The general formulation of this method is explained by using Bregman distance. The Bregman distance associated with a convex function 𝐽 at the point 𝑣 is 𝐷^𝑝_𝐽(𝑢,𝑣)=𝐽(𝑢)−𝐽(𝑣)−𝑝,𝑢−𝑣, (6) where 𝑝 is in the subgradient of 𝐽 at 𝑣. Clearly, this is not a distance in the usual sense, it measures closeness in the sense that 𝐷^𝑝_𝐽(𝑢,𝑣)≥0 and 𝐷^𝑝_𝐽(𝑢,𝑣)≥𝐷^𝑝_𝐽(𝑤,𝑣) for 𝑤 on the line segment between 𝑢 and 𝑣. Again, consider two convex energy functionals, 𝐽 and 𝐿, defined over 𝑅^𝑛 with min_𝑢∈𝑅^𝑛𝐿(𝑢)=0. The associated unconstrained minimization problem is min_𝑢𝐽(𝑢)+𝜇𝐿(𝑢). (7) Modify this problem by iteratively solving 𝑢^𝑘+1=min_𝑢𝐷^𝑝_𝐽𝑢,𝑢^𝑘(𝑢),𝑢+𝜇𝐿_𝑘+1=min_𝑢𝐽(𝑢)−𝑝^𝑘,𝑢−𝑢^𝑘+𝜇𝐿(𝑢) (8)(9) as suggested by Bregman in [21]. For simplicity, assume that 𝐿 is differentiable, and in this case, we have 0∈𝛿(𝐷^𝑝_𝐽(𝑢,𝑢^𝑘)+𝜇𝐿(𝑢)), where this subdifferential is evaluated at 𝑢^𝑘+1. Since 𝑝^𝑘+1∈𝛿𝐽(𝑢^𝑘+1) at this location, then 𝑝^𝑘+1=𝑝^𝑘−∇𝐿𝑢^𝑘+1. (10) Apply the Bregman iteration (8) on (5) 𝑢^𝑘+1=min_𝑢𝐷^𝑝_𝐽𝑢,𝑢^𝑘+𝜇2‖𝐾𝑢−𝑣‖²₂,𝑢^𝑘+1=min_𝑢𝐽(𝑢)−𝑝^𝑘,𝑢−𝑢^𝑘+𝜇2‖𝐾𝑢−𝑣‖²₂,𝑝^𝑘+1=𝑝^𝑘−𝜇𝐾^𝑇𝐾𝑢^𝑘+1−𝑣. (11) Bregman iterations of this form were considered in [12]. Here, it is shown that when 𝐾 is linear, this seemingly complicated iteration is equivalent to the simplified method 𝑢^𝑘+1=min_𝑢𝐽(𝑢)+𝜇2𝐾𝑢−𝑣^𝑘2₂,𝑣^𝑘+1=𝑣^𝑘+𝑣−𝐾𝑢^𝑘+1. (12)

2.2. Split Bregman Method

Apply Bregman framework to solve 𝐿₁-regularized problem (4). In this case, assume that both Φ(𝑢) and 𝐿(𝑢) are convex functions and Φ(𝑢) is differentiable. In this method, decouple the 𝐿₁ and 𝐿₂ parts of energy in (4). Now, consider (4) as min_𝑢,𝑑|𝑑|+𝐿(𝑢)suchthat𝑑=Φ(𝑢). (13) To solve this problem, first convert it into an unconstrained problem min_𝑢,𝑑|𝑑|(𝑢)+𝜇+𝐿2‖𝑑−Φ(𝑢)‖²₂. (14) Let 𝐽(𝑢,𝑑)=|𝑑|+𝐿(𝑢), and let 𝐾(𝑢,𝑑)=𝑑−Φ(𝑢), then we can see that (14) is simply an application of (5). Use (11) to obtain the solution for the above problem: 𝑢^𝑘+1,𝑑^𝑘+1=min_𝑢,𝑑𝐷^𝑝_𝐽𝑢,𝑢^𝑘,𝑑,𝑑^𝑘+𝜇2‖𝑑−Φ(𝑢)‖²₂,𝑢^𝑘+1,𝑑^𝑘+1=min_𝑢,𝑑𝐽(𝑢,𝑑)−𝑝^𝑘_𝑢,𝑢−𝑢^𝑘−𝑝^𝑘_𝑑,𝑑−𝑑^𝑘+𝜇2‖𝑑−Φ(𝑢)‖²₂,𝑝_𝑢^𝑘+1=𝑝^𝑘_𝑢−𝜇(∇Φ)^𝑇Φ𝑢^𝑘+1−𝑑^𝑘+1,𝑝_𝑑^𝑘+1=𝑝^𝑘_𝑑−𝜇𝑑^𝑘+1−Φ𝑢^𝑘+1. (15) Applying simplification presented in (12), we get the following solutions: 𝑢^𝑘+1,𝑑^𝑘+1=min_𝑢,𝑑|𝑑|(𝑢)+𝜇+𝐿2𝑑−Φ(𝑢)−𝑣^𝑘2₂,𝑣^𝑘+1=𝑣^𝑘+Φ𝑢^𝑘+1−𝑑^𝑘+1. (16)(17) In order to perform the minimization effectively, split the 𝐿₁ and 𝐿₂ components of (16) and minimize with respect to 𝑢 and 𝑑 separately. The two steps result in the following solutions: 𝑢^𝑘+1=min_𝑢𝐿(𝑢)+𝜇2𝑑^𝑘−Φ(𝑢)−𝑣^𝑘2₂,𝑑^𝑘+1=min_𝑑|𝑑|+𝜇2𝑑−Φ𝑢^𝑘+1−𝑣^𝑘2₂. (18)(19) In (19), there is no coupling between elements of 𝑑. We can use shrinkage operator to find the optimal value of 𝑑 as follows: 𝑑^𝑘+1=shrinkΦ(𝑢)_𝑗+𝑣^𝑘_𝑗,1𝜇, (20) where ()=𝑥shrink𝑥,𝜆|𝑥|∗max(|𝑥|−𝜆,0). (21) This shrinkage operation is extremely fast and requires only few operations per element of 𝑑^𝑘+1.

Based on the above described equations, the generalized Split Bregman algorithm can be explained in Algorithm 1.

while  𝑢𝑘−𝑢𝑘−12>𝑡𝑜𝑙  do  for  𝑛=1 to 𝑁  do   𝑢𝑘+1=min𝑢𝐿(𝑢)+𝜇2𝑑𝑘−Φ(𝑢)−𝑣𝑘22   𝑑𝑘+1=min𝑑|𝑑|+𝜇2𝑑−Φ𝑢𝑘+1−𝑣𝑘22  end for  𝑣𝑘+1=𝑣𝑘+Φ𝑢𝑘+1−𝑑𝑘+1 end while

Algorithm 1 

Generalized Split Bregman algorithm.

2.3. Nonlocal Total Variation Norm

Total Variation (TV) regularization [24, 25] makes the recovered image quality sharper, but they do not preserve the fine structures, details, and textures. This effect is caused by the regularity assumption of the TV formulation of the image model, namely, that the image has a simple geometric description consisting of a set of connected sets (objects) with smooth contours (edges). Additionally, the model assumes that the image is smooth inside single objects and has discontinuous jumps across the boundaries. Therefore, TV regularization is optimal to reduce the noise and to reconstruct the main geometrical configuration in an image. However, it fails to preserve texture, details, and fine structures, because they behave in all aspects like noise and thus cannot be distinguished from noise. The nonlocal total variation (NLTV) extends the TV functional to a nonlocal variant using the definition of nonlocal derivative operators based on a nonlocal weight function [14–17, 26, 27]. NLTV is an effective tool instead of TV for improving the signal-to-noise ratio in practical application [14–16]. Recently, it has been successfully used for 4D computed tomography reconstruction from few projection data [28]. NLTV extends the TV functional to a nonlocal variant using the definition of nonlocal derivative operators based on a nonlocal weight function (graph). The notion nonlocal means that any point can directly interact with any other point in the whole image domain, where the intensity of the interaction is depending on the value of the weight function. This weight function should represent the similarity of the two points and should be significant, if both points are similar in an appropriate measure. Therefore, the expectation is that such an approach is able to process both structures (geometrical parts) and texture within the same framework, due to the identification of recurring structures in the whole image. A brief review regarding the definition of nonlocal functional is given below.

Let Ω∈𝑅², and let 𝑢∶Ω→𝑅 be a real function. Assume that 𝑤∶Ω×Ω→𝑅 is a weight function which is symmetric and nonnegative. Then, the nonlocal gradient of an image 𝑢(𝑥) is defined as ∇_NL𝑢𝑥,𝑦∶=𝑢(𝑥)−𝑢𝑦𝑤𝑥,𝑦∶Ω×Ω⟶𝑅. (22) The norm of a vector 𝑞∶Ω×Ω→𝑅 at point 𝑥∈Ω is given by 𝑞(𝑥)=_Ω𝑞𝑥,𝑦²𝑑𝑦. (23) Hence, the norm of the nonlocal gradient of a function 𝑢∶Ω→𝑅 at point 𝑥∈Ω is represented as ∇_NL𝑢(𝑥)=_Ω𝑢(𝑥)−𝑢𝑦²𝑤𝑥,𝑦𝑑𝑦∶Ω⟶𝑅₊. (24) The nonlocal divergence operator can be defined using the standard adjoint relation with the nonlocal gradient as follows: ∇_NL𝑢,𝑞∶=−𝑢,div_NL𝑞,∀𝑢∶Ω⟶𝑅,𝑞∶Ω×Ω⟶𝑅 (25) which guides to the definition of nonlocal divergence of the nonlocal vector 𝑞div_NL𝑞(𝑥)=_Ω𝑞𝑥,𝑦−𝑞𝑦,𝑥𝑤𝑥,𝑦𝑑𝑦∶Ω⟶𝑅. (26) Next, the nonlocal Laplacian operator is defined as Δ_NL𝑢(𝑥)1∶=2div_NL∇_NL𝑢(𝑥=_Ω)𝑢(𝑥)−𝑢𝑦𝑤𝑥,𝑦𝑑𝑦. (27) The discrete forms of the nonlocal gradient, divergence, and laplace operators can be represented as follows: ∇_NL𝑢_𝑖∶=_𝑗𝑢_𝑗−𝑢_𝑖²𝑤_𝑖𝑗,div_NL𝑞_𝑖∶=_𝑗𝑞_𝑖,𝑗−𝑞_𝑗,𝑖𝑤_𝑖𝑗,Δ_NL𝑢_𝑖∶=_𝑗𝑢_𝑗−𝑢_𝑖𝑤_𝑖𝑗, (28) where 𝑢_𝑖 represents the value of a pixel 𝑖 in the image (1≤𝑖≤𝑁) and 𝑤_𝑖𝑗 is the discrete sparse version of 𝑤(𝑥,𝑦), and it is defined as 𝑤−𝐺𝑥,𝑦=exp^𝜎∗𝑢(𝑥+⋅)−𝑢𝑦+⋅²2ℎ², (29) where ℎ is a filtering parameter; in general ℎ corresponds to the noise level; usually, we set it to be the standard deviation of the noise, and 𝐺_𝜎 is a Gaussian of standard deviation 𝜎, and 𝑢(𝑥) and 𝑢(𝑦) are the image values in pixel 𝑥 and 𝑦. The weight functions 𝑤(𝑥,𝑦) denote how much the difference between pixels 𝑥 and 𝑦 is penalized in the images. The more similar the neighborhoods of 𝑥 and 𝑦 are, more the difference should be penalized, and vice versa. The weight functions 𝑤(𝑥,𝑦), significant only if the window around 𝑦 looks like the corresponding window around 𝑥. Hence, the nonlocal algorithm is very efficient in reducing noise, while preserving contrast in natural images and redundant structures such as texture. In our work, we used 5×5 pixel patches, a search neighborhood window of size 11×11. The observed noisy data 𝑣 is taken as the reference image to construct the weight, and by this weighed averaging, the structures, for example, boundaries, are reinforced, while the noise gets cancelled out [29]. The weights are computed by using either a distance between the noisy pixel values |𝑢(𝑥)−𝑢(𝑦)| [30–32] or a distance between the patches around 𝑥 and 𝑦 [33–35]. The use of NLTV reduces the noise in the reconstructed image, thus the difference between reference and reconstructed image reduces. From the definition of signal-to-noise ratio (SNR), it is clear that the reduction in image difference increases the SNR.

3. Proposed Work

The proposed work jointly minimizes a linear combination of nonlocal total variation and least-square data-fitting term to reconstruct the MR image from undersampled data. The main aim is to solve the compressed sensing MRI problem (2) using Split Bregman algorithm and nonlocal total variation. In this work, the nonlocal total variation is taken as the 𝐿₁-regularization functional and solved using Split Bregman iteration. Recall (2) min_𝑢𝐽(𝑢)suchthat‖𝐾𝑢−𝑣‖²₂<𝜎². (30) Using Bregman iteration method, (30) can be reduced to a sequence of unconstrained problems of the form 𝑢^𝑘+1=min_𝑢𝐽(𝑢)+𝜇2𝐾𝑢−𝑣^𝑘2₂,𝑣^𝑘+1=𝑣^𝑘+𝑣−𝐾𝑢^𝑘+1, (31)(32) where 𝐽(𝑢) represents 𝐿₁-regularization term. In order to proceed further selection of regularization method is important. Here, we choose nonlocal total variation (NLTV) as the regularizer; that is 𝐽(𝑢)=∇_NL𝑢. (33) Now, (31) becomes min_𝑢∇_NL𝑢₁+𝜇2𝐾𝑢−𝑣^𝑘2₂, (34) where ‖∇_NL𝑢‖₁=∑_𝑖|∇_NL𝑢|_𝑖. We can write (34) as follows by introducing an auxiliary variable 𝑑 instead of ∇_NL𝑢: min_𝑢,𝑑‖𝑑‖₁+𝜇2𝐾𝑢−𝑣^𝑘2₂suchthat𝑑=∇_NL𝑢. (35) Equation (35) can be converted into unconstrained form by using the quadratic penalty method min_𝑢,𝑑‖𝑑‖₁+𝛾2𝑑−∇_NL𝑢²₂+𝜇2𝐾𝑢−𝑣^𝑘2₂. (36) Using split Bregman method, (36) can be transformed into the following forms: min_𝑢,𝑑‖𝑑‖₁+𝛾2𝑑−∇_NL𝑢−𝑏^𝑘2₂+𝜇2𝐾𝑢−𝑣^𝑘2₂,𝑏^𝑘+1=𝑏^𝑘+∇_NL𝑢^𝑘+1−𝑑^𝑘+1. (37)(38) Equation (37) is convex and can be minimized by alternatively solving the following two minimization subproblems with respect to 𝑢 and 𝑑𝑢^𝑘+1=argmin_𝑢𝛾2𝑑^𝑘−∇_NL𝑢−𝑏^𝑘2₂+𝜇2𝐾𝑢−𝑣^𝑘2₂,𝑑^𝑘+1=argmin_𝑑‖𝑑‖₁+𝛾2𝑑−∇_NL𝑢−𝑏^𝑘2₂. (39) By direct computation, the optimal conditions of (39) are −𝛾div_NL𝑑^𝑘−∇_NL𝑢−𝑏^𝑘+𝜇𝐾^𝑇𝐾𝑢−𝑣^𝑘𝑑=0,_𝑘+1=∇_NL𝑢^𝑘+1+𝑏^𝑘∇_NL𝑢^𝑘+1+𝑏^𝑘max∇_NL𝑢^𝑘+1+𝑏^𝑘−1𝛾,0. (40)(41) Use the Gauss-Seidel iteration to get a fast solution of (40), and the discrete solution is represented as 𝑢_𝑖^𝑘+1=1𝐾^𝑇𝐾𝜇+𝛾∑_𝑗𝑤_𝑖𝑗×𝛾_𝑗𝑤_𝑖𝑗𝑢^𝑘+𝜇𝑣^𝑘_𝑖+𝛾_𝑗𝑤_𝑖𝑗𝑑^𝑘_𝑖𝑗−𝑏^𝑘_𝑖𝑗−𝑑^𝑘_𝑗𝑖+𝑏^𝑘_𝑗𝑖. (42) As explained in the introduction, 𝐾 is partial Fourier matrix (𝐾=𝑃ℱ), where 𝑃∈𝑅^𝑀×𝑁 is an identity matrix (𝑀≪𝑁) and ℱ is a discrete Fourier matrix. Using the identity ℱ^𝑇=ℱ⁻¹, now (42) becomes 𝑢_𝑖^𝑘+1=1ℱ⁻¹𝑃^𝑇𝑃ℱ𝜇+𝛾∑_𝑗𝑤_𝑖𝑗×𝛾_𝑗𝑤_𝑖𝑗𝑢^𝑘+𝜇𝑣^𝑘_𝑖+𝛾_𝑗𝑤_𝑖𝑗𝑑^𝑘_𝑖𝑗−𝑏^𝑘_𝑖𝑗−𝑑^𝑘_𝑗𝑖+𝑏^𝑘_𝑗𝑖. (43) The discrete solution for (41) can be written as 𝑑^𝑘+1_𝑖,𝑗=𝑤_𝑖𝑗𝑢_𝑗^𝑘+1−𝑢_𝑖^𝑘+1+𝑏^𝑘_𝑖𝑗∑_𝑗𝑤_𝑖𝑗𝑢_𝑗^𝑘+1−𝑢_𝑖^𝑘+12+𝑏^𝑘_𝑖𝑗²×max_𝑗𝑤_𝑖𝑗𝑢_𝑗^𝑘+1−𝑢_𝑖^𝑘+12+𝑏^𝑘_𝑖𝑗²−1𝛾,0. (44)

Finally, the Bregman variable is updated as 𝑏^𝑘+1_𝑖𝑗=𝑏^𝑘_𝑖𝑗+_𝑗𝑤_𝑖𝑗𝑢_𝑗^𝑘+1−𝑢_𝑖^𝑘+1−𝑑^𝑘+1_𝑖𝑗. (45) The proposed method is summarized as Algorithm 2.

Initialize: 𝑢0=ℱ−1𝑣,𝑑0=𝑏0=0,𝑘=0 and 𝜇,𝛾>0 while  𝑢𝑘−𝑢𝑘−1>𝑡𝑜𝑙  do  Compute 𝑢𝑘+1 based on (43)  Compute 𝑑𝑘+1 based on (44)  Update 𝑏𝑘 based on (45)  𝑘=𝑘+1 end while

Algorithm 2 

Proposed algorithm.

4. Evaluation of Image Quality

In this work, a detailed evaluation study has done on the reconstruction of MR images, which represent varying degrees of object structural complexity. Even though algorithms based on regularization techniques effectively remove streaks, other aspects of image quality should also be analyzed. To address this, a number of image quality evaluations are performed at different levels including qualitative visualization-based evaluation and quantitative metric-based evaluation.

4.1. Qualitative Visualization-Based Evaluation

In qualitative visualization-based evaluation, reconstructed image obtained with different algorithms are visually compared with the reference image.

4.2. Quantitative Metric-Based Evaluation

Besides the visualization-based evaluation, similarity between reconstructed and reference images is quantitatively assessed by means of four measures such as signal-to-noise ratio (SNR), relative error (RE), structural similarity index (SSIM), and feature similarity index (FSIM). SNR and RE are widely used for measuring reconstruction accuracy, SSIM and FSIM are used for evaluating the similarity between reconstructed and reference image.

4.2.1. Signal-to-Noise Ratio (SNR)

One can see that 𝑢SNR=10log^ref−𝑢𝑢_ref−𝑢_rec, (46) where 𝑢_ref is the reference image, 𝑢 is the mean intensity value of 𝑢_ref, and 𝑢_rec is the reconstructed image.

4.2.2. Relative Error (RE)

One can see that 𝑢RE=^rec−𝑢_ref𝑢_ref×100%. (47)

4.2.3. Structural Similarity Index (SSIM)

The SNR measurement gives a numerical value on the damage, but it does not describe its type. Moreover, as is noted in [36, 37], it does not quite represent the quality perceived by human observers. For medical imaging applications, where images are degraded, must eventually be examined by experts, traditional evaluation remains insufficient. For this reason, objective approaches are needed to assess the medical imaging quality. We then evaluate a new paradigm to estimate the quality of medical images based on the assumption that the human visual system (HVS) is highly adapted to extract structural information. The similarity index compares the brightness 𝐼(𝑥,𝑦), contrast 𝑐(𝑥,𝑦), and structure 𝑠(𝑥,𝑦) between each pair of vectors, where the SSIM index between two signals 𝑥 and 𝑦 is given by the following expression [38, 39]: SSIM𝑥,𝑦=𝐼𝑥,𝑦𝑐𝑥,𝑦𝑠𝑥,𝑦. (48) However, the comparison of brightness is determined by the following expression: 𝐼𝑥,𝑦=2𝜇_𝑥𝜇_𝑦+𝐶₁𝜇_𝑥+𝜇_𝑦+𝐶₁, (49) where the average intensity of signal 𝑥 is given by 𝜇_𝑥=1𝑁^𝑁_𝑖=1𝑥_𝑖,𝐶₁=𝐾₁𝐿² (50) the constant 𝐾₁≪1, and 𝐿 is the dynamic row of the pixel values (255 for an image coded on 8 bits). The function of contrast comparison takes the following form: 𝑐𝑥,𝑦=2𝜎_𝑥𝜎_𝑦𝜎²_𝑥+𝜎²_𝑦+𝐶₂, (51) where 𝜎_𝑥=𝜇_𝑥(𝑥²)−𝜇²_𝑥(𝑥) is the standard deviation of the original signal 𝑥, 𝐶₂=(𝐾₂𝐿)², and the constant 𝐾₂≪1.

The function of structure comparison is defined as follows: 𝑠=𝜎𝑥,𝑦^𝑥𝑦+𝐶₃𝜎_𝑥𝜎_𝑦+𝐶₃=cov𝑥,𝑦+𝐶₃𝜎_𝑥𝜎_𝑦+𝐶₃, (52) where cov(𝑥,𝑦)=𝜇_𝑥𝑦−𝜇_𝑥𝜇_𝑦 and 𝐶₃=𝐶₂/2.

Then, the expression of the structural similarity index becomes SSIM𝑥,𝑦=2𝜇_𝑥𝜇_𝑦+𝐶₁2𝜎_𝑥𝑦+𝐶₂𝜇²_𝑥+𝜇²_𝑦+𝐶₁𝜎²_𝑥+𝜎²_𝑦+𝐶₂. (53)

4.2.4. Feature Similarity Index (FSIM)

SSIM index provides image quality assessment (IQA) from pixel-based stage to structure-based stage. Human visual system (HVS) understands an image mainly based on its low-level features: mainly, the phase congruency (PC), which is a measure of the significance of a local structure and it is dimensionless. PC is used as the primary feature in FSIM [40]. The secondary feature used in FSIM is the image gradient magnitude (GM). In order to find out the feature similarity between two images 𝑓₁ and 𝑓₂, the above mentioned parameters PC and GM are to be calaculated first.

Let PC₁, PC₂, 𝐺₁, and 𝐺₂ be the phase congruency and gradient magnitude of images 𝑓₁ and 𝑓₂, respectively. Initially, separate the feature similarity measurement between 𝑓₁(𝑥) and 𝑓₂(𝑥) into two components. First similarity measure is based on PC₁(𝑥) and PC₂(𝑥) and is defined as 𝑆_PC(𝑥)=2PC₁(𝑥)⋅PC₂(𝑥)+𝑇₁PC²₁(𝑥)+PC²₂(𝑥)+𝑇₁, (54) where 𝑇₁ is a positive constant which increases the stability of 𝑆_PC. Value of 𝑇₁ depends on the dynamic range of PC. Similarly, the similarity measure based on GM values 𝐺₁(𝑥) and 𝐺₂(𝑥) is defined as 𝑆_𝐺(𝑥)=2𝐺₁(𝑥)⋅𝐺₂(𝑥)+𝑇₂𝐺²₁(𝑥)+𝐺²₂(𝑥)+𝑇₂, (55) where 𝑇₁ is a positive constant which depends on the dynamic range of GM value.

Next step is to combine 𝑆_PC(𝑥) and 𝑆_𝐺(𝑥) to get the similarity measure 𝑆_𝐿(𝑥) between 𝑓₁(𝑥) and 𝑓₂(𝑥) and is defined as 𝑆_𝐿(𝑥)=𝑆_PC(𝑥)^𝛼⋅𝑆_𝐺(𝑥)^𝛽. (56) The relative importance of PC and GM features can be adjusted by means of the parameters 𝛼 and 𝛽. For simplicity, set 𝛼=𝛽=1, then (56) becomes 𝑆_𝐿(𝑥)=𝑆_PC(𝑥)⋅𝑆_𝐺(𝑥), (57) where 𝑆_𝐿(𝑥) represents the similarity at each location 𝑥, and the overall similarity should be found. For a given location 𝑥, if any of 𝑓₁(𝑥) and 𝑓₂(𝑥) has a significant PC value, it implies that this position x will have a high impact on HVS in evaluating the similarity between 𝑓₁ and 𝑓₂. Therefore, introducing a new term PC_𝑚(𝑥)=max(PC₁(𝑥),PC₂(𝑥)) to weight the importance of 𝑆_𝐿(𝑥) in the overall similarity between 𝑓₁ and 𝑓₂. Accordingly, the FSIM index can defined as ∑FSIM=^𝑥∈Ω𝑆_𝐿(𝑥)⋅PC_𝑚(𝑥)∑_𝑥∈ΩPC_𝑚(𝑥), (58) where Ω means the whole image spatial domain.

5. Experiments and Numerical Results

The experimental setup used in previous works [5–7] is explained here. Suppose that an MR image 𝑢 has 𝑚×𝑚 pixels and the partial Fourier transform 𝐾 in (3) consists of 𝑛 rows of a 𝑚×𝑚 matrix corresponding to the full 2D discrete Fourier transform. The 𝑛 selected rows correspond to the acquired data 𝑣. The sampling ratio is outlined as 𝑛/𝑚. The scanning time is shorter if the sampling ratio is smaller. In the 𝑘-space, randomly choose more samples in low frequencies and fewer samples in higher frequencies. This sample theme has been widely used for compressed MR image reconstruction, and therefore similar themes are utilized in [4–7]. Practically, the sampling scheme and speed in MR imaging also depend on the physical and physiological limitations [4]. In the proposed work, the compressive sensing matrix 𝐾=𝑃ℱ, where 𝑃 is a row selector matrix, and ℱ is the Fourier transform matrix. For an 𝑁×𝑁 image, we randomly choose m coefficients, then 𝑅 is a sampling matrix of size 𝑚×𝑁². All experiments are conducted on a PC with an Intel core-i72670, 2.2 GHz CPU in MATLAB environment. The result of the proposed work is compared with the existing methods like TVCMRI [5], RecPF [6], CSA, FCSA [7], and SB-TV [8]. The observation measurement 𝑣 is synthesized as 𝑣=𝐾𝑢+𝑛, where 𝑛 is the noise, 𝐾 and 𝑣 are given as inputs, and 𝑢 is the unknown target. In this work we considered two sets of observations: one is contaminated with Gaussian noise of standard deviation 𝜎=0.01, and the other is contaminated with Rician noise of noise level 5%.

The proposed and existing algorithms are tested using four 2D MR images: brain, chest, artery, and cardiac, imges respectively, as shown in Figure 1. They have been used for validation in [7]. For convenience, they have the same size of 256×256. The sample ratio is set to be approximately 20%. To perform comparisons, all methods run 50 iterations, except SB-TV and the proposed method. SB-TV completes reconstruction in 35 iterations, while the proposed method takes only 30 iterations. The proposed method provides best visual effects on all MR images. Figures 2, 3, 4, and 5 show the visual comparisons of the reconstructed results using Gaussian noisy observations by different methods in the brain, chest, artery, and cardiac images, respectively. Figure 10 gives the performance comparisons between different methods in terms of the CPU time over SNR. Table 1 shows the average value of CPU time, SNR, RE, SSIM and FSIM of different methods in Gaussian noise case. Figures 6, 7, 8, and 9 show the visual comparisons of the reconstructed results using Rician noisy observations by different methods in the brain, chest, artery, and cardiac images, respectively. Figure 11 gives the performance comparisons between different methods in terms of the CPU time over SNR. Table 2 shows the average value of CPU time, SNR, RE, SSIM and FSIM of different methods in Rician noise case. In both cases, the reconstructed results produced by the proposed method is better than those produced by the TVCMRI, RecPF, CSA, FCSA, and SB-TV. The reconstruction performance of the proposed work is the best in terms of both the computational complexity and reconstruction accuracy, which clearly demonstrate the efficiency of this method for the compressed MR image reconstruction.

Table 1 

Comparison of the average values of quality evaluation parameters (Gaussian noise case).

Table 2 

Comparison of the average values of quality evaluation parameters (Rician noise case).

(a)  Brain

(b) Chest

(c) Artery

(d) Cardiac

Figure 1 

2D MR test images.

(a) TVCMRI (SNR=14.19 dB)

(b) ecPF (SNR=14.75 dB)

(c) CSA (SNR=15.20 dB)

(d) FCSA (SNR=15.78 dB)

(e) SB-TV (SNR=16.31)

(f) Proposed (SNR=17.49 dB)

Figure 2 

Brain MR image reconstruction from 20% samples (Gaussian noise case).

(a) TVCMRI (SNR=16.99 dB)

(b) RecPF (SNR=17.38 dB)

(c) CSA (SNR=17.13 dB)

(d) FCSA (SNR=17.63 dB)

(e) SB-TV (SNR=18.03)

(f) Proposed (SNR=18.89 dB)

Figure 3 

Chest MR image reconstruction from 20% samples (Gaussian noise case).

(a) TVCMRI (SNR=18.70 dB)

(b) RecPF (SNR=19.98 dB)

(c) CSA (SNR=22.15 dB)

(d) FCSA (SNR=23.61 dB)

(e) SB-TV (SNR=24.17)

(f) Proposed (SNR=25.75 dB)

Figure 4 

Artery MR image reconstruction from 20% samples (Gaussian noise case).

(a) TVCMRI (SNR=16.97 dB)

(b) RecPF (SNR=17.65 dB)

(c) CSA (SNR=17.95 dB)

(d) FCSA (SNR=18.77 dB)

(e) SB-TV (SNR=18.98)

(f) Proposed (SNR=19.47 dB)

Figure 5 

Cardiac MR image reconstruction from 20% samples (Gaussian noise case).

(a) TVCMRI (SNR=13.31 dB)

(b) RecPF (SNR=13.85 dB)

(c) CSA (SNR=14.81 dB)

(d) FCSA (SNR=15.31 dB)

(e) SB-TV (SNR=16.24)

(f) Proposed (SNR=17.08dB)

Figure 6 

Brain MR image reconstruction from 20% samples (Rician noise case).

(a) TVCMRI (SNR=15.30 dB)

(b) RecPF (SNR=15.50 dB)

(c) CSA (SNR=15.45dB)

(d) FCSA (SNR=15.69 dB)

(e) SB-TV (SNR=16.67)

(f) Proposed (SNR=17.18 dB)

Figure 7 

Chest MR image reconstruction from 20% samples (Rician noise case).

(a) TVCMRI (SNR=18.34 dB)

(b) RecPF (SNR=19.55 dB)

(c) CSA (SNR=21.26 dB)

(d) FCSA (SNR=22.54 dB)

(e) SB-TV (SNR=23.64)

(f) Proposed (SNR=24.28 dB)

Figure 8 

Artery MR image reconstruction from 20% samples (Rician noise case).

(a) TVCMRI (SNR=14.68 dB)

(b) RecPF (SNR=14.95 dB)

(c) CSA (SNR=14.98 dB)

(d) FCSA (SNR=15.26dB)

(e) SB-TV (SNR=17.52)

(f) Proposed (SNR=18.93 dB)

Figure 9 

Cardiac MR image reconstruction from 20% samples (Rician noise case).

(a) Brain image

(b) Chest image

(c) Artery image

(d) Cardiac image

Figure 10 

Performance comparisons (CPU time versus SNR) on different MR images (Gaussian noise case).

(a) Brain image

(b) Chest image

(c) Artery image

(d) Cardiac image

Figure 11 

Performance comparisons (CPU time versus SNR) on different MR images (Rician noise case).

6. Conclusion

An efficient algorithm is proposed for the compressed MR image reconstruction based on nonlocal total variation and Split Bregman method. The algorithm effectively solves a NLTV-based 𝐿₁-norm term by using the Split Bregman method. NLTV can effectively avoid blocky artifacts caused by traditional TV regularization. Numerous experiments were conducted to compare different reconstruction methods. The results of our study indicate that the proposed method outperforms the classic methods in terms of both accuracy and complexity.

Acknowledgment

The authors would like to thank the Government of Canada for the financial support through the Commonwealth Scholarship 2011-12 for this research work.