Computational and Mathematical Methods in Medicine

Volume 2013 (2013), Article ID 985819, 16 pages

http://dx.doi.org/10.1155/2013/985819

Research Article

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

Varun P. Gopi,¹ P. Palanisamy,¹ Khan A. Wahid,² and Paul Babyn³

¹Department of Electronics and Communication Engineering, National Institute of Technology, Tiruchirappalli, Tamil Nadu 620015, India

²Department of Electrical and Computer Engineering, University of Saskatchewan, Saskatoon, SK, Canada

³Department of Medical Imaging, Royal University Hospital, University of Saskatchewan, Saskatoon, SK, Canada

Received 26 April 2013; Revised 26 June 2013; Accepted 11 July 2013

Academic Editor: Dimitri Van De Ville

Copyright © 2013 Varun P. Gopi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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

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

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

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

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

where is in the subgradient of at . Clearly, this is not a distance in the usual sense, it measures closeness in the sense that and for on the line segment between and . Again, consider two convex energy functionals, and , defined over with . The associated unconstrained minimization problem is

Modify this problem by iteratively solving

as suggested by Bregman in [21]. For simplicity, assume that is differentiable, and in this case, we have , where this subdifferential is evaluated at . Since at this location, then

Apply the Bregman iteration (8) on (5)

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

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

To solve this problem, first convert it into an unconstrained problem

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:

Applying simplification presented in (12), we get the following solutions:

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:

In (19), there is no coupling between elements of . We can use shrinkage operator to find the optimal value of as follows:

where

This shrinkage operation is extremely fast and requires only few operations per element of .

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

while    do

    for   to   do

    end for

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

The norm of a vector at point is given by

Hence, the norm of the nonlocal gradient of a function at point is represented as

The nonlocal divergence operator can be defined using the standard adjoint relation with the nonlocal gradient as follows:

which guides to the definition of nonlocal divergence of the nonlocal vector

Next, the nonlocal Laplacian operator is defined as

The discrete forms of the nonlocal gradient, divergence, and laplace operators can be represented as follows:

where represents the value of a pixel in the image () and is the discrete sparse version of , and it is defined as

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 pixel patches, a search neighborhood window of size . 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)

Using Bregman iteration method, (30) can be reduced to a sequence of unconstrained problems of the form

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

Now, (31) becomes

where . We can write (34) as follows by introducing an auxiliary variable instead of :

Equation (35) can be converted into unconstrained form by using the quadratic penalty method

Using split Bregman method, (36) can be transformed into the following forms:

Equation (37) is convex and can be minimized by alternatively solving the following two minimization subproblems with respect to and

By direct computation, the optimal conditions of (39) are

Use the Gauss-Seidel iteration to get a fast solution of (40), and the discrete solution is represented as

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

The discrete solution for (41) can be written as

Finally, the Bregman variable is updated as

The proposed method is summarized as Algorithm 2.

Initialize: and

while    do

    Compute based on (43)

    Compute based on (44)

    Update based on (45)

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

where is the reference image, is the mean intensity value of , and is the reconstructed image.

4.2.2. Relative Error (RE)

One can see that

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]:

However, the comparison of brightness is determined by the following expression:

where the average intensity of signal is given by

the constant , 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:

where is the standard deviation of the original signal , , and the constant .

The function of structure comparison is defined as follows:

where and .

Then, the expression of the structural similarity index becomes

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 , , , 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 and and is defined as

where is a positive constant which increases the stability of . Value of depends on the dynamic range of PC. Similarly, the similarity measure based on GM values and is defined as

where is a positive constant which depends on the dynamic range of GM value.

Next step is to combine and to get the similarity measure between and and is defined as

The relative importance of PC and GM features can be adjusted by means of the parameters and . For simplicity, set , then (56) becomes

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 to weight the importance of in the overall similarity between and . Accordingly, the FSIM index can defined as

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 , 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 . 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).

Image
TVCMRI
RecPF
CSA
FCSA
SB-TV
Proposed

SNR (dB)
Brain
14.19
14.75
15.20
15.78
16.31
17.49

Chest
16.99
17.38
17.13
17.63
18.03
18.89

Artery
18.70
19.98
22.15
23.61
24.17
25.75

Cardiac
16.97
17.65
17.95
18.77
18.98
19.47

RE ()
Brain
16.96
16.54
07.97
07.37
06.90
05.98

Chest
11.58
11.21
10.84
10.33
08.89
07.96

Artery
12.61
11.84
05.77
04.90
04.57
04.10

Cardiac
14.03
13.39
11.10
09.78
08.29
07.72

SSIM
Brain
0.9567
0.9742
0.9906
0.9935
0.9940
0.9947

Chest
0.9829
0.9839
0.9848
0.9861
0.9887
0.9914

Artery
0.9825
0.9835
0.9848
0.9919
0.9924
0.9940

Cardiac
0.9614
0.9642
0.9748
0.9799
0.9835
0.9889

FSIM
Brain
0.8712
0.8773
0.9489
0.9654
0.9719
0.9765

Chest
0.8990
0.9012
0.9063
0.9148
0.9237
0.9374

Artery
0.9010
0.9168
0.9516
0.9668
0.9705
0.9742

Cardiac
0.9313
0.9366
0.9454
0.9551
0.9587
0.9629

CPU time (s)
Brain
2.14
2.05
1.96
1.88
1.78
1.62

Chest
2.21
2.16
1.95
1.86
1.73
1.62

Artery
2.26
2.15
1.86
1.83
1.72
1.64

Cardiac
2.13
2.10
1.84
1.82
1.80
1.67

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

Image
TVCMRI
RecPF
CSA
FCSA
SB-TV
Proposed

SNR (dB)
Brain
13.31
13.85
14.81
15.31
16.24
17.08

Chest
15.30
15.50
15.45
15.69
16.67
17.18

Artery
18.34
19.55
21.26
22.54
23.64
24.28

Cardiac
14.68
14.95
14.98
15.26
17.52
18.93

RE ()
Brain
17.15
16.78
13.01
12.39
10.36
10.01

Chest
13.19
12.94
11.58
10.98
09.38
08.81

Artery
13.24
12.46
06.76
05.95
05.04
04.71

Cardiac
15.06
14.73
12.69
10.24
09.89
08.96

SSIM
Brain
0.9511
0.9625
0.9787
0.9807
0.9860
0.9870

Chest
0.9747
0.9762
0.9801
0.9816
0.9828
0.9889

Artery
0.9810
0.9824
0.9832
0.9889
0.9904
0.9924

Cardiac
0.9582
0.9592
0.9697
0.9717
0.9812
0.9858

FSIM
Brain
0.8499
0.8562
0.8994
0.9082
0.9290
0.9339

Chest
0.8860
0.8968
0.9031
0.9067
0.9221
0.9298

Artery
0.8652
0.8797
0.9424
0.9566
0.9672
0.9715

Cardiac
0.9220
0.9304
0.9352
0.9404
0.9518
0.9565

CPU time (s)
Brain
2.10
2.05
1.96
1.88
1.78
1.62

Chest
2.26
2.14
1.97
1.85
1.73
1.63

Artery
2.18
2.14
1.88
1.84
1.74
1.63

Cardiac
2.10
2.12
1.84
1.82
1.79
1.67

(a)  Brain

(b) Chest

(c) Artery

(d) Cardiac

Figure 1: 2D MR test images.

(a) TVCMRI ( dB)

(b) ecPF ( dB)

(c) CSA ( dB)

(d) FCSA ( dB)

(e) SB-TV ()

(f) Proposed ( dB)

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

(a) TVCMRI ( dB)

(b) RecPF ( dB)

(c) CSA ( dB)

(d) FCSA ( dB)

(e) SB-TV ()

(f) Proposed ( dB)

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

(a) TVCMRI ( dB)

(b) RecPF ( dB)

(c) CSA ( dB)

(d) FCSA ( dB)

(e) SB-TV ()

(f) Proposed ( dB)

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

(a) TVCMRI ( dB)

(b) RecPF ( dB)

(c) CSA ( dB)

(d) FCSA ( dB)

(e) SB-TV ()

(f) Proposed ( dB)

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

(a) TVCMRI ( dB)

(b) RecPF ( dB)

(c) CSA ( dB)

(d) FCSA ( dB)

(e) SB-TV ()

(f) Proposed (dB)

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

(a) TVCMRI ( dB)

(b) RecPF ( dB)

(c) CSA (dB)

(d) FCSA ( dB)

(e) SB-TV ()

(f) Proposed ( dB)

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

(a) TVCMRI ( dB)

(b) RecPF ( dB)

(c) CSA ( dB)

(d) FCSA ( dB)

(e) SB-TV ()

(f) Proposed ( dB)

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

(a) TVCMRI ( dB)

(b) RecPF ( dB)

(c) CSA ( dB)

(d) FCSA (dB)

(e) SB-TV ()

(f) Proposed ( 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.

References

E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489–509, 2006.

D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.

L. He, T. C. Chang, S. Osher, T. Fang, and P. Speier, “MR image reconstruction by using the iterative refinement method and nonlinear inverse scale space methods,” UCLA CAM Report, 2006.

M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: the application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182–1195, 2007.

S. Ma, W. Yin, Y. Zhang, and A. Chakraborty, “An efficient algorithm for compressed MR imaging using total variation and wavelets,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR '08), pp. 1–8, Anchorage, Alaska, USA, June 2008.

J. Yang, Y. Zhang, and W. Yin, “A fast alternating direction method for TVL1-L2 signal reconstruction from partial Fourier data,” IEEE Journal on Selected Topics in Signal Processing, vol. 4, no. 2, pp. 288–297, 2010.

J. Huang, S. Zhang, and D. Metaxas, “Efficient MR image reconstruction for compressed MR imaging,” Medical Image Analysis, vol. 15, no. 5, pp. 670–679, 2011.

T. Goldstein and S. Osher, “The split Bregman method for L1-regularized problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323–343, 2009.

D. C. Dobson and F. Santosa, “Recovery of blocky images from noisy and blurred data,” SIAM Journal on Applied Mathematics, vol. 56, no. 4, pp. 1181–1198, 1996.

A. Chambolle and P. L. Lions, “Image recovery via total variation minimization and related problems,” Numerische Mathematik, vol. 76, no. 2, pp. 167–188, 1997.

Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution Equations: The Fifteenth Dean Jacqueline B. Lewis Memorial Lectures, vol. 22 of University Lecture Series, American Mathematical Society, Boston, Mass, USA, 2001.

S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, “An iterative regularization method for total variation-based image restoration,” Multiscale Modeling and Simulation, vol. 4, no. 2, pp. 460–489, 2005.

Z. Tian, X. Jia, K. Yuan, T. Pan, and S. B. Jiang, “Low-dose CT reconstruction via edge-preserving total variation regularization,” Physics in Medicine and Biology, vol. 56, no. 18, pp. 5949–5967, 2011.

G. Gilboa and S. Osher, “Nonlocal operators with applications to image processing,” Multiscale Modeling and Simulation, vol. 7, no. 3, pp. 1005–1028, 2008.

X. Zhang, M. Burger, X. Bresson, and S. Osher, “Bregmanized nonlocal regularization for deconvolution and sparse reconstruction,” SIAM Journal on Imaging Sciences, vol. 3, no. 3, pp. 253–276, 2010.

G. Peyre, S. Bougleux, and L. Cohen, “Non-local regularization of inverse problems,” in Computer Vision—ECCV 2008, vol. 5304 of Lecture Notes in Computer Science, pp. 57–68, 2008.

Y. Lou, X. Zhang, S. Osher, and A. Bertozzi, “Image recovery via nonlocal operators,” Journal of Scientific Computing, vol. 42, no. 2, pp. 185–197, 2010.

X. Zhang, M. Burger, X. Bresson, and S. Osher, “Bregmanized nonlocal regularization for deconvolution and sparse reconstruction,” SIAM Journal on Imaging Sciences, vol. 3, no. 3, pp. 253–276, 2010.

S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, Cambridge, UK, 2004.

J. Nocedal and S. Wright, Numerical Optimization, Springer, New York, NY, USA, 2nd edition, 2006.

L. M. Bregman, “The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming,” USSR Computational Mathematics and Mathematical Physics, vol. 7, no. 3, pp. 200–217, 1967.

A. Osher, Y. Mao, B. Dong, and W. Yin, “Fast linearized Bregman iterations for compressed sensing and sparse denoising,” UCLA CAM Report, University of California at Los Angeles, Los Angeles, Calif, USA, 2008.

L. He, T. C. Chang, and S. Osher, “MR image reconstruction from sparse radial samples by using iterative refinement procedures,” in Proceedings of the 13th Annual Meeting of ISMRM, p. 696, 2005.

L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D, vol. 60, no. 1–4, pp. 259–268, 1992.

A. Chambolle, “An algorithm for total variation minimization and applications,” Journal of Mathematical Imaging and Vision, vol. 20, no. 1-2, pp. 89–97, 2004.

M. Mignotte, “A non-local regularization strategy for image deconvolution,” Pattern Recognition Letters, vol. 29, no. 16, pp. 2206–2212, 2008.

A. Elmoataz, O. Lezoray, and S. Bougleux, “Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing,” IEEE Transactions on Image Processing, vol. 17, no. 7, pp. 1047–1060, 2008.

X. Jia, Y. Lou, B. Dong, Z. Tian, and S. Jiang, “4D computed tomography reconstruction from few-projection data via temporal non-local regularization,” in Medical Image Computing and Computer-Assisted Intervention—MICCAI 2010, vol. 6361 of Lecture Notes in Computer Science, pp. 143–150, 2010.

A. Buades, B. Coll, and J. M. Morel, “A review of image denoising algorithms, with a new one,” Multiscale Modeling and Simulation, vol. 4, no. 2, pp. 490–530, 2005.

A. Spira, R. Kimmel, and N. Sochen, “A short-time Beltrami kernel for smoothing images and manifolds,” IEEE Transactions on Image Processing, vol. 16, no. 6, pp. 1628–1636, 2007.

C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” in Proceedings of the 6th IEEE International Conference on Computer Vision, pp. 839–846, Bombay, India, January 1998.

L. P. Yaroslavsky, Digital Picture Processing—An Introduction, Springer, Berlin, Germany, 1985.

R. R. Coifman, S. Lafon, A. B. Lee et al., “Geometric diffusions as a tool for harmonic analysis and structure definition of data: diffusion maps,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 21, pp. 7426–7431, 2005.

A. Buades, B. Coll, and J. M. Morel, “A review of image denoising algorithms, with a new one,” Multiscale Modeling and Simulation, vol. 4, no. 2, pp. 490–530, 2005.

A. Szlam, M. Maggioni, and R. R. Coifman, “A general framework for adaptive regularization based on diffusion processes on graphs,” Journal of Machine Learning Research, vol. 9, pp. 1711–1739, 2008.

W. S. Geisler and M. S. Banks, “Visual performance,” in Handbook of Optics, vol. 1, McGraw-Hill, New York, NY, USA, 1995.

A. B. Watson and L. B. Kreslake, “Measurement of visual impairment scales for digital video,” in 6th Human Vision and Electronic Imaging, vol. 4299 of Proceedings of SPIE, pp. 79–89, San Jose, Calif, USA, January 2001.

Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004.

Z. Wang and A. C. Bovik, “A universal image quality index,” IEEE Signal Processing Letters, vol. 9, no. 3, pp. 81–84, 2002.

L. Zhang, L. Zhang, X. Mou, and D. Zhang, “FSIM: a feature similarity index for image quality assessment,” IEEE Transactions on Image Processing, vol. 20, no. 8, pp. 2378–2386, 2011.