# Median Filter Based Compressed Sensing Model with Application to MR Image Reconstruction

Authors: Yunyun Yang; Xuxu Qin; Boying Wu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Magnetic resonance imaging (MRI) has become a helpful technique and developed rapidly in clinical medicine and diagnosis. Magnetic resonance (MR) images can display more clearly soft tissue structures and are important for doctors to diagnose diseases. However, the long acquisition and transformation time of MR images may limit their application in clinical diagnosis. Compressed sensing methods have been widely used in faithfully reconstructing MR images and greatly shorten the scanning and transforming time. In this paper we present a compressed sensing model based on median filter for MR image reconstruction. By combining a total variation term, a median filter term, and a data fitting term together, we first propose a minimization problem for image reconstruction. The median filter term makes our method eliminate additional noise from the reconstruction process and obtain much clearer reconstruction results. One key point of the proposed method lies in the fact that both the total variation term and the median filter term are presented in the L1 norm formulation. We then apply the split Bregman technique for fast minimization and give an efficient algorithm. Finally, we apply our method to numbers of MR images and compare it with a related method. Reconstruction results and comparisons demonstrate the accuracy and efficiency of the proposed model.


---

## Body

## 1. Introduction

Compressed sensing [1–3] is a new sampling theory appearing in signal and image processing communities; it allows us to directly acquire a data with compressed representation. In other words, the compressed sensing method can reconstruct the original signal accurately by using a small number of linear combinations of the compression transform coefficients. As a promising method, the theory has been introduced to image reconstruction [4–7], wireless sensor networks [8–10], speech coding [11, 12], image memorability prediction [13], and others [14–16]. Its broad scope and generality has made the CS method become a hot research topic.In clinical medicine and diagnosis, MRI technology has developed quickly and become a useful technique. MR images have the advantage that soft tissue structures can be displayed more clearly. However, due to the defects of imaging equipment and technology, the acquisition time of MR images may be long. In MRI, compressed sensing has been used to shorten magnetic resonance imaging scanning sessions on conventional hardware.One classic CS model based on the total variation norm was proposed in [1], which is mainly used for piecewise constant images. But this model could not reconstruct piecewise smooth images well, another CS model based on the total variation norm, and the Besov norm was proposed in [5] to overcome this difficulty. By adding a wavelet transform, this model improved the reconstruction accuracy for inhomogeneous images, however, it is much more time-consuming. The split Bregman technique [17] is an efficient method which can be successfully used in image denoising, segmentation and compressed sensing [18–21]. The authors in [21] applied this technique into the above CS model [5] and gave a split Bregman algorithm to improve the efficiency with promising results. In this paper, we call the CS model based on the split Bregman algorithm as the SBA model. The SBA model can obtain relatively satisfactory reconstruction results to some extent.Besides the long acquisition time, MR images often contain noises, which may be from the acquisition or transmission process and degrade the quality of MR images. Therefore, image filtering or image smoothing is also a necessary step in MR image pre-processing. Many mature image filtering methods come out, such as Gabor filter, Gaussian filter, Wiener filter, and the median filter. The median filter is an effective nonlinear smoothing method, which can remove isolated noises and preserve details of the image boundary.Therefore, in this paper we combine the median filter into the traditional compressed sensing reconstruction technique to propose a new CS model. The proposed model can reconstruct MR images with high precision and eliminate isolated noise from reconstruction procedure. Furthermore, the proposed energy functional has the special structure of theL1 norm; thus the split Bregman technique can be used for fast minimization. Then we test our model with many MR images and compare the reconstruction results with those of the SBA model qualitatively and quantitatively. Experimental results validate the superiority of our model in terms of accuracy and efficiency.The rest of the paper is organized as follows. In Section2 we briefly review related concepts and methods. In Section 3 we give the main contribution of this paper. We first present the proposed minimization problem of our method, and then we apply the split Bregman technique to efficiently minimize the minimization problem. In Section 4, first we give a fast algorithm, and then we apply our method to test numbers of brain MR images and present experimental results. We conclude this paper in Section 5.


## 2. Related Concepts and Methods

### 2.1. Concept of Compressed Sensing

The notion of compressed sensing as put forward in the works of Donoho [2] and others [1, 3] has suggested that we can accurately reconstruct a sparse signal x from fewer measurements than the traditional method needs; we mean that the amount of sampled data by the compressed sensing method is much smaller than the number of sampled data specified by the Nyquist sampling law, but the signal also can be reconstructed very accurately. In some applications, the CS helps to reduce the sample rate and save the system resources. The samples are linear incoherent and can be described as follows:(1)y=Ax,where A∈Rm×n is the measurement matrix, usually m≪n, y∈Rm is the measurement vector, and x∈Rn is a sparse signal. Note that the measurement matrix A is usually chosen as the random matrices [22, 23].The reconstructedx can be realized by solving the following L0-minimization problem:(2)minx⁡x0s.t.y=Ax,where the definition of x0 is the number of nonzero elements in vector x. However, Natarajan [24] proved that the problem (2) is a nondeterministic polynomial hard (NP-hard) problem in general. As an alternative, the convex relaxation of (2) has provided an approximation solution to x [25, 26], which leads to the following convex optimization problem:(3)minx⁡x1s.t.y=Ax.Since the observed signaly usually contains noise, thereby the Lagrangian version of the problem (3) would be more suitable:(4)minx⁡x1s.t.y-Ax22<θ,where θ is expected to the variance of noise.The aboveL1 minimization problem is a nonsmooth optimization problem; thus many efficient optimization methods are also available, such as [21, 27, 28]. In the real world, many kinds of signals are not sparse but they usually have a sparse representation under a suitable basis, for example, the wavelet transform or the Fourier transform. It is known that the k-space of MR images is encoded with the Fourier transform; thus MR images are naturally compressible in the frequency domain. The MR image restoration problems from the partial k-space data will be introduced in Section 2.2.


### 2.2. The Model Based on Total Variation Norm

Candès et al. [1] presented the MR image restoration problem under the total variation (TV) regularization as(5)minv⁡▿v1+σ2RFv-Sp22,where Sp represents the observed partial k-space data, R represents the measurement matrix, and σ is a positive parameter. F is the Fourier transform and Fv is used to simulate the full k-space data. The total variation, also known as the TV norm, is constructed as follows:(6)▿v1=∑i,j▿vx,i,j2+▿vy,i,j2.This model is easy to understand and implement, but it has a disadvantage that it is mainly for piecewise constant images. It has been shown that models based on only the TV norm may have low accuracy in reconstructing inhomogeneous images [4, 5].


### 2.3. The Model Based on Total Variation Norm and Besov Norm

Most MR images are considered as piecewise smooth; thus in the image domain the TV norm can be used for sparse representation. Actually, it has been claimed that in the frequency domain the underlying images can also be sparse [5]. The authors applied the Haar orthogonal wavelet transform Ψ and combined the Besov norm into the following minimization problem:(7)minv⁡Ψv1+α▿v1+σ2RFv-Sp22,where α is a positive parameter.This model combined the total variation norm and the Besov norm together; thus it has higher reconstruction accuracy for piecewise smooth images. But the reconstruction efficiency might be decreased because adding the Besov norm with wavelet transform is more time-consuming.


## 2.1. Concept of Compressed Sensing

The notion of compressed sensing as put forward in the works of Donoho [2] and others [1, 3] has suggested that we can accurately reconstruct a sparse signal x from fewer measurements than the traditional method needs; we mean that the amount of sampled data by the compressed sensing method is much smaller than the number of sampled data specified by the Nyquist sampling law, but the signal also can be reconstructed very accurately. In some applications, the CS helps to reduce the sample rate and save the system resources. The samples are linear incoherent and can be described as follows:(1)y=Ax,where A∈Rm×n is the measurement matrix, usually m≪n, y∈Rm is the measurement vector, and x∈Rn is a sparse signal. Note that the measurement matrix A is usually chosen as the random matrices [22, 23].The reconstructedx can be realized by solving the following L0-minimization problem:(2)minx⁡x0s.t.y=Ax,where the definition of x0 is the number of nonzero elements in vector x. However, Natarajan [24] proved that the problem (2) is a nondeterministic polynomial hard (NP-hard) problem in general. As an alternative, the convex relaxation of (2) has provided an approximation solution to x [25, 26], which leads to the following convex optimization problem:(3)minx⁡x1s.t.y=Ax.Since the observed signaly usually contains noise, thereby the Lagrangian version of the problem (3) would be more suitable:(4)minx⁡x1s.t.y-Ax22<θ,where θ is expected to the variance of noise.The aboveL1 minimization problem is a nonsmooth optimization problem; thus many efficient optimization methods are also available, such as [21, 27, 28]. In the real world, many kinds of signals are not sparse but they usually have a sparse representation under a suitable basis, for example, the wavelet transform or the Fourier transform. It is known that the k-space of MR images is encoded with the Fourier transform; thus MR images are naturally compressible in the frequency domain. The MR image restoration problems from the partial k-space data will be introduced in Section 2.2.


## 2.2. The Model Based on Total Variation Norm

Candès et al. [1] presented the MR image restoration problem under the total variation (TV) regularization as(5)minv⁡▿v1+σ2RFv-Sp22,where Sp represents the observed partial k-space data, R represents the measurement matrix, and σ is a positive parameter. F is the Fourier transform and Fv is used to simulate the full k-space data. The total variation, also known as the TV norm, is constructed as follows:(6)▿v1=∑i,j▿vx,i,j2+▿vy,i,j2.This model is easy to understand and implement, but it has a disadvantage that it is mainly for piecewise constant images. It has been shown that models based on only the TV norm may have low accuracy in reconstructing inhomogeneous images [4, 5].


## 2.3. The Model Based on Total Variation Norm and Besov Norm

Most MR images are considered as piecewise smooth; thus in the image domain the TV norm can be used for sparse representation. Actually, it has been claimed that in the frequency domain the underlying images can also be sparse [5]. The authors applied the Haar orthogonal wavelet transform Ψ and combined the Besov norm into the following minimization problem:(7)minv⁡Ψv1+α▿v1+σ2RFv-Sp22,where α is a positive parameter.This model combined the total variation norm and the Besov norm together; thus it has higher reconstruction accuracy for piecewise smooth images. But the reconstruction efficiency might be decreased because adding the Besov norm with wavelet transform is more time-consuming.


## 3. The Proposed Model and Fast Minimization

### 3.1. The Proposed Minimization Problem

In this section, we present the following minimization problem:(8)minv▿v1s.t.RFv-Sp22<δ,Mvk-v1<μ.where δ and μ are two positive parameters taking small values. M is a two-dimensional median filter. Given an original image v(x,y), the output image u(x,y) is defined as(9)ux,y=Mvx,y=medianvx-k,y-l,k,l,W,where W is a template whose size is usually 3×3 or 5×5.The proposed minimization problem (8) is then converted into an unconstrained problem as follows:(10)minv⁡▿v1+Mvk-v1+σ2RFv-Sp22,where σ>0 is a parameter.We explain the proposed minimization problem of our model in detail. From the minimization problem (5) of the model based on the TV norm, we know that v-vk→0 when the algorithm converges, where vk is the last iteration result. In our model, vk is replaced with Mvk and another TV norm Mvk-v1 is added into the problem (5). Mvk is the result by applying the median filtering M to the kth iteration reconstruction image, which can reduce the additional noise that comes from the reconstruction process, and Mvk-v1 will force v to approach Mvk in each iteration, which means we can also achieve the effect of eliminating noise in the image v. Therefore, by jointly minimizing the TV term ▿v1, the data fitting term RFv-Sp22, and the median filter term Mvk-v1, we can obtain a much clearer reconstructed image.The median filter termMvk-v1 in the proposed model works by alleviating the aliasing produced in the iterative process. With this term, the proposed model shows more precise reconstructed results than the SBA model. Moreover, the median filter term can also reduce the cumulative error of the algorithm and accelerate the convergence of the algorithm.


### 3.2. Fast Minimization with the Split Bregman Method

Note that the proposed minimization problem (10) is indeed a L1 minimization problem; thus we can apply the split Bregman method to solve the problem (10). We first introduce two auxiliary variables d→ and g→ and obtain the following minimization problem:(11)minv,d→,g→⁡d→1+g→1+σ2RFv-Sp22+η2d→-▿v22+β2g→-Mvk-v22,where η and β are two positive parameters. The two quadratic functions (η/2)d→-▿v22 and (β/2)g→-(Mvk-v)22 are used to enforce weakly the constraints d→=▿v and g→=Mvk-v. In order to establish strong constraints to the above two equations, we need to introduce additional relaxation variables b→ and bg→ andthen the problem becomes the following sequence problems:(12)minv,d→,g→,b→,bg→⁡d→1+g→1+σ2RFv-Sp22+η2d→-▿v-b→k+β2g→-Mvk-v-bg→k22,where we update the Bregman variables b→ and bg→ as follows:(13)b→k+1=b→k+▿vk+1-d→k+1,bg→k+1=bg→k+Mvk-vk+1-g→k+1.We apply the alternating minimization scheme to solve the problem (12) with respect to v, d→ and g→ as follows:(14)vk+1=argminv⁡σ2RFv-Sp22+η2d→k-▿v-b→k+β2g→k-Mvk-v-bg→k22,(15)d→k+1=argmind→⁡d→1+η2d→-▿vk+1-b→k22,(16)g→k+1=argming→⁡g→1+β2g→-Mvk+vk+1-bg→k22.For the minimization problem (14), we can obtain the following equation with respect to v:(17)σFTRTRF-η▵+βIvk+1=sk,where(18)sk=σFTRSp+η▿·d→k-b→k+βg→k-Mvk-bg→k.Since the Fourier transform is orthogonal, we have FT=F-1, let K=(uRTR-ηFT▵F+βI), and then get(19)vk+1=F-1K-1Fsk.For the other minimization problems (15) and (16), we can use the shrinkage operator to explicitly give the updating schemes for the variables d→ and g→ as follows:(20)d→k+1=shrink▿vk+1+b→k,1η,(21)g→k+1=shrinkMvk-vk+1+bg→k,1β,where the vector shrinkage operator is defined as(22)shrinkx→,ζ=x→x→max⁡x→-ζ,0→,x→≠0→,0→,x→=0→.


## 3.1. The Proposed Minimization Problem

In this section, we present the following minimization problem:(8)minv▿v1s.t.RFv-Sp22<δ,Mvk-v1<μ.where δ and μ are two positive parameters taking small values. M is a two-dimensional median filter. Given an original image v(x,y), the output image u(x,y) is defined as(9)ux,y=Mvx,y=medianvx-k,y-l,k,l,W,where W is a template whose size is usually 3×3 or 5×5.The proposed minimization problem (8) is then converted into an unconstrained problem as follows:(10)minv⁡▿v1+Mvk-v1+σ2RFv-Sp22,where σ>0 is a parameter.We explain the proposed minimization problem of our model in detail. From the minimization problem (5) of the model based on the TV norm, we know that v-vk→0 when the algorithm converges, where vk is the last iteration result. In our model, vk is replaced with Mvk and another TV norm Mvk-v1 is added into the problem (5). Mvk is the result by applying the median filtering M to the kth iteration reconstruction image, which can reduce the additional noise that comes from the reconstruction process, and Mvk-v1 will force v to approach Mvk in each iteration, which means we can also achieve the effect of eliminating noise in the image v. Therefore, by jointly minimizing the TV term ▿v1, the data fitting term RFv-Sp22, and the median filter term Mvk-v1, we can obtain a much clearer reconstructed image.The median filter termMvk-v1 in the proposed model works by alleviating the aliasing produced in the iterative process. With this term, the proposed model shows more precise reconstructed results than the SBA model. Moreover, the median filter term can also reduce the cumulative error of the algorithm and accelerate the convergence of the algorithm.


## 3.2. Fast Minimization with the Split Bregman Method

Note that the proposed minimization problem (10) is indeed a L1 minimization problem; thus we can apply the split Bregman method to solve the problem (10). We first introduce two auxiliary variables d→ and g→ and obtain the following minimization problem:(11)minv,d→,g→⁡d→1+g→1+σ2RFv-Sp22+η2d→-▿v22+β2g→-Mvk-v22,where η and β are two positive parameters. The two quadratic functions (η/2)d→-▿v22 and (β/2)g→-(Mvk-v)22 are used to enforce weakly the constraints d→=▿v and g→=Mvk-v. In order to establish strong constraints to the above two equations, we need to introduce additional relaxation variables b→ and bg→ andthen the problem becomes the following sequence problems:(12)minv,d→,g→,b→,bg→⁡d→1+g→1+σ2RFv-Sp22+η2d→-▿v-b→k+β2g→-Mvk-v-bg→k22,where we update the Bregman variables b→ and bg→ as follows:(13)b→k+1=b→k+▿vk+1-d→k+1,bg→k+1=bg→k+Mvk-vk+1-g→k+1.We apply the alternating minimization scheme to solve the problem (12) with respect to v, d→ and g→ as follows:(14)vk+1=argminv⁡σ2RFv-Sp22+η2d→k-▿v-b→k+β2g→k-Mvk-v-bg→k22,(15)d→k+1=argmind→⁡d→1+η2d→-▿vk+1-b→k22,(16)g→k+1=argming→⁡g→1+β2g→-Mvk+vk+1-bg→k22.For the minimization problem (14), we can obtain the following equation with respect to v:(17)σFTRTRF-η▵+βIvk+1=sk,where(18)sk=σFTRSp+η▿·d→k-b→k+βg→k-Mvk-bg→k.Since the Fourier transform is orthogonal, we have FT=F-1, let K=(uRTR-ηFT▵F+βI), and then get(19)vk+1=F-1K-1Fsk.For the other minimization problems (15) and (16), we can use the shrinkage operator to explicitly give the updating schemes for the variables d→ and g→ as follows:(20)d→k+1=shrink▿vk+1+b→k,1η,(21)g→k+1=shrinkMvk-vk+1+bg→k,1β,where the vector shrinkage operator is defined as(22)shrinkx→,ζ=x→x→max⁡x→-ζ,0→,x→≠0→,0→,x→=0→.


## 4. Numerical Algorithm and Experimental Results

### 4.1. Numerical Algorithm

Based on the above minimization schemes presented in Section3.1, we give the following fast algorithm.Algorithm 1.
Step 1. Input R and Sp, initialize v0=0, Mv0=0, d→0=g→0=0, and b→0=bg→0=0, and give values for η, σ, β, and μ.
while v k + 1 - v k > ρ do
Step 2. Update vk+1 by (19).
Step 3. Update d→k+1 by (20).
Step 4. Update g→k+1 by (21).
Step 5. Update b→k+1,bg→k+1 by (13).
end while
Step 6. Output the reconstruct result vk+1.In Algorithm1, the parameter ρ is the tolerance criteria. The parameters η=0.5, σ=0.5, β=0.005, and ρ=0.001 are chosen for the algorithm. It is worth noting that values of these parameters are empirical values summarized in our experiments. According to our experiences, the parameter η should be set the same as the parameter σ, and we follow the rule that the higher the resolution of the image is, the greater the values of the parameters η and σ should be taken. The parameter β is used to control the median filter term. Generally, when the value of the parameter β is less than 0.01, our model can obtain accurate reconstruction results. If the value of ρ is too large, the precision of reconstruction will be reduced. If the value of ρ is too small, for example, 0.0001, the reconstruction time will be greatly increased, but the accuracy of reconstruction will not be significantly improved. Therefore, a more moderate tolerance is recommended in our algorithm.


### 4.2. Experimental Results

In this section we apply our model to lots of MR images and present experimental results. We also provide the comparison results between our model and the SBA model to demonstrate the improvement of our model. Both algorithms are run on a Lenove desktop with Intel(R) Core(TM)i5-6500 CPU, 3.20GHz 4GB RAM, with Matlab 2014a on Windows 10 operating system.In terms of the reconstruction accuracy, an evaluation value called the peak signal-to-noise ratio (PSNR) is used in this paper. If we denote the mean square error between the original image and the processed image by MSE, then the PSNR value is defined as(23)PSNR=20×log10⁡255MSE.From this definition, larger PSNR value implies higher reconstruction accuracy. The relative error (Rerror) is another index to evaluate the accuracy of reconstructed images, which is defined as(24)Rerror=vk-v22v22,where vk is the reconstruction image and v is the reference image.In Figure1 we present the sampling mask used in our experiments. This kind of radial sampling mask is designed to simulate the sampling method used in medical MR imaging. The under sampling ratio is defined as(25)Samp.Ratio=sumRnumelR,where sum(R) is the total number of the selected sampling pixels, while nembel(R) is the total number of pixels in R.Figure 1
The radial sampling mask R. The size of R is equal to that of the underlying images, where white pixels represent the positions of sampling pixels.In Figure2 we give the reconstruction results of a synthetic 2D phantom image with different sampling ratios, where (a) shows the reference image and (b), (c), and (d) are the recovered images wiht the sampling ratios of 3.9%, 5.8%, and 9.5%, respectively. This figure shows that our algorithm can give an acceptable reconstruction result of the phantom image with a very low sampling ratio of 3.9% and the sampling ratio of 9.5% is sufficient for reconstructing a faithful image.Figure 2
2D phantom image reconstruction with different sampling ratios. (a) The reference image; (b), (c), and (d) the recovered images with the sampling ratios of 3.9%, 5.8%, and 9.5%, respectively.
(a) (b) (c) (d)In Figure3 we present the reconstruction results of a real 2D brain image with different sampling ratios, where the source of this image can be found in [29]. We denote it by Brain 1 in Tables 1 and 2. Figure 3(a) shows the reference image, and we can see that the structure of this brain image is not very complex; hence 35.01% of the sampled data is enough to reconstruct an image with satisfactory accuracy.Table 1
Numerical results for the synthetic phantom image and real MR images with different sampling ratios.
Image Samp.Ratio PSNR (dB) Time (s) Rerror Phantom 3.90% 28.19 6.34 0.165 512×512 5.76% 34.55 4.65 0.078 9.50% 42.49 3.08 0.040  Image Samp.Ratio PSNR (dB) Time (s) Rerror  Brain 1 7.60% 29.89 1.74 0.120 512×512 15.00% 34.04 1.55 0.070 35.01% 37.47 1.55 0.051  Image Samp.Ratio PSNR (dB) Time (s) Rerror  Brain 2 7.60% 31.05 1.59 0.130 512×512 18.50% 37.52 1.55 0.060 35.00% 44.46 1.53 0.020Table 2
Numerical comparison for different MR images between our model and the SBA model.
Image Method Samp.Ratio PSNR (dB) Time (s) Rerror Brain 1 SBA 25.35% 34.63 5.32 0.067 512×512 Our model 25.35% 36.06 1.58 0.053  Brain 2 SBA 35.00% 40.04 8.47 0.050 512×512 Our model 35.00% 43.84 1.57 0.032  Brain 3 SBA 35.00% 36.34 1.44 0.029 200×200 Our model 35.00% 41.67 0.19 0.015Figure 3
2D brain image reconstruction with different sampling ratios. (a) The reference image; (b), (c), and (d) the recovered images with the sampling ratios of 7.6%, 15.0%, and 35.1%, respectively.
(a) (b) (c) (d)Figure4 shows that our model can also recover a faithful image with the sampling ratio of 35.0% of a real brain MR image with a higher level of complexity. The reference brain MR image is from the paper [30], and we denote it by Brain 2 in Tables 1 and 2.Figure 4
2D brain image reconstruction with different sampling ratios. (a) The reference image; (b), (c), and (d) the recovered images with the sampling ratios of 7.6%, 18.5%, and 35.0%, respectively.
(a) (b) (c) (d)In Figure5 we present the comparison results of our method and the SBA method in reconstructing three different real brain MR images, in which the image (i) can be found in [31] and denoted by Brain 3 in this paper. Column 1 shows the sample data, where image (a) is with the sampling ratio of 25.35% and images (e) and (i) are both with sampling ratios of 35.01%. Column 2 and Column 3 show the reconstruction results of the proposed model and the SBA model, respectively. The last column shows the reference images. We can see from Figure 5 that some fine structures of the reference image are missing with the SBA model, while our model can reconstruct clearer structures, which can be obviously observed by comparing (f) and (g).Figure 5
The comparison results of our method and the SBA method for three different brain MR images. Column 1: The sample data. Column 2: The reconstruction results of the proposed model. Column 3: The reconstruction results of the SBA model. Column 4: The reference brain MR images.In Figure6 we show the comparison results of our method and the SBA method with different sampling ratios in terms of the PSNR value and the relative error, where the results of the phantom image, the image Brain 1 and the image Brain 2 are shown in Column 1, Column 2, and Column 3, respectively. The first row of Figure 6 shows the PSNR values while the second row shows the relative errors. From Figure 6 we can see that PSNR values of the recovered images with our model are higher than those with the SBA model, and the relative errors are much smaller than those of the SBA model. This indicates that our model performs better in recovering MR images from undersampled images than the SBA model.Figure 6
The comparison results of our method and the SBA method with different sampling ratios. Row 1: The comparison of PSNR (dB) with the proposed model and the SBA model. Row 2: The comparison of relative errors with the proposed model and the SBA model.We evaluate the improvement of the proposed model by comparing the PSNR value and the computational time with the SBA model for a number of 100 brain MR images. The comparison results are shown in Figure7. Note that the sources of these brain images are the same as Brain 3 and the sizes of all images are 200×200. From Figure 7(a), we can observe that the PSNR values of reconstructed images with our model are all greater than those with the SBA model, which indicates the higher reconstruction accuracy of our model. Besides, from Figure 7(b), it can be seen the average computational time with our model is below 0.2s, much less than that with the SBA model, which is above 1.4s. We can conclude that our model is much more efficient and accurate than the SBA model.Figure 7
Comparison of reconstruction results between the proposed model and the SBA model with application to 100 brain MR images. (a) The PSNR (dB); (b) the computational time (s).
(a) (b)


### 4.3. Quantitative Analysis

In this section, to clearly see the numerical comparison of the reconstruction accuracy and efficiency with different sampling ratios for different images, we record the size, PSNR value, computational time, and relative error for the three images in Figures2, 3, and 4 in Table 1. From Table 1 we can observe that for the same image, with a higher sampling ratio, the PSNR value is larger and the Rerror value is smaller, which means that the reconstruction result is more accurate, while the CPU time is relatively less and thus is more efficient. This an intuitive and understandable result.We compare the numerical comparison results of our method and the SBA method for the three real brain MR images in Table2. Column 4 shows that the PSNR values of our model are 3.52 dB higher than those of the SBA model on average, which validates the improvement of reconstruction accuracy. Column 5 shows the computation time of the SBA model and the proposed model, from which we can obtain that the computation time of our model is much less than that of the SBA model. Besides, the Rerror values shown in Column 6 also verify higher reconstruction precision of our model than the other one. Table 2 demonstrates that our model can obtain reconstruction results with higher accuracy but less time. Therefore, our model is more accurate and efficient compared with the SBA method.


## 4.1. Numerical Algorithm

Based on the above minimization schemes presented in Section3.1, we give the following fast algorithm.Algorithm 1.
Step 1. Input R and Sp, initialize v0=0, Mv0=0, d→0=g→0=0, and b→0=bg→0=0, and give values for η, σ, β, and μ.
while v k + 1 - v k > ρ do
Step 2. Update vk+1 by (19).
Step 3. Update d→k+1 by (20).
Step 4. Update g→k+1 by (21).
Step 5. Update b→k+1,bg→k+1 by (13).
end while
Step 6. Output the reconstruct result vk+1.In Algorithm1, the parameter ρ is the tolerance criteria. The parameters η=0.5, σ=0.5, β=0.005, and ρ=0.001 are chosen for the algorithm. It is worth noting that values of these parameters are empirical values summarized in our experiments. According to our experiences, the parameter η should be set the same as the parameter σ, and we follow the rule that the higher the resolution of the image is, the greater the values of the parameters η and σ should be taken. The parameter β is used to control the median filter term. Generally, when the value of the parameter β is less than 0.01, our model can obtain accurate reconstruction results. If the value of ρ is too large, the precision of reconstruction will be reduced. If the value of ρ is too small, for example, 0.0001, the reconstruction time will be greatly increased, but the accuracy of reconstruction will not be significantly improved. Therefore, a more moderate tolerance is recommended in our algorithm.


## 4.2. Experimental Results

In this section we apply our model to lots of MR images and present experimental results. We also provide the comparison results between our model and the SBA model to demonstrate the improvement of our model. Both algorithms are run on a Lenove desktop with Intel(R) Core(TM)i5-6500 CPU, 3.20GHz 4GB RAM, with Matlab 2014a on Windows 10 operating system.In terms of the reconstruction accuracy, an evaluation value called the peak signal-to-noise ratio (PSNR) is used in this paper. If we denote the mean square error between the original image and the processed image by MSE, then the PSNR value is defined as(23)PSNR=20×log10⁡255MSE.From this definition, larger PSNR value implies higher reconstruction accuracy. The relative error (Rerror) is another index to evaluate the accuracy of reconstructed images, which is defined as(24)Rerror=vk-v22v22,where vk is the reconstruction image and v is the reference image.In Figure1 we present the sampling mask used in our experiments. This kind of radial sampling mask is designed to simulate the sampling method used in medical MR imaging. The under sampling ratio is defined as(25)Samp.Ratio=sumRnumelR,where sum(R) is the total number of the selected sampling pixels, while nembel(R) is the total number of pixels in R.Figure 1
The radial sampling mask R. The size of R is equal to that of the underlying images, where white pixels represent the positions of sampling pixels.In Figure2 we give the reconstruction results of a synthetic 2D phantom image with different sampling ratios, where (a) shows the reference image and (b), (c), and (d) are the recovered images wiht the sampling ratios of 3.9%, 5.8%, and 9.5%, respectively. This figure shows that our algorithm can give an acceptable reconstruction result of the phantom image with a very low sampling ratio of 3.9% and the sampling ratio of 9.5% is sufficient for reconstructing a faithful image.Figure 2
2D phantom image reconstruction with different sampling ratios. (a) The reference image; (b), (c), and (d) the recovered images with the sampling ratios of 3.9%, 5.8%, and 9.5%, respectively.
(a) (b) (c) (d)In Figure3 we present the reconstruction results of a real 2D brain image with different sampling ratios, where the source of this image can be found in [29]. We denote it by Brain 1 in Tables 1 and 2. Figure 3(a) shows the reference image, and we can see that the structure of this brain image is not very complex; hence 35.01% of the sampled data is enough to reconstruct an image with satisfactory accuracy.Table 1
Numerical results for the synthetic phantom image and real MR images with different sampling ratios.
Image Samp.Ratio PSNR (dB) Time (s) Rerror Phantom 3.90% 28.19 6.34 0.165 512×512 5.76% 34.55 4.65 0.078 9.50% 42.49 3.08 0.040  Image Samp.Ratio PSNR (dB) Time (s) Rerror  Brain 1 7.60% 29.89 1.74 0.120 512×512 15.00% 34.04 1.55 0.070 35.01% 37.47 1.55 0.051  Image Samp.Ratio PSNR (dB) Time (s) Rerror  Brain 2 7.60% 31.05 1.59 0.130 512×512 18.50% 37.52 1.55 0.060 35.00% 44.46 1.53 0.020Table 2
Numerical comparison for different MR images between our model and the SBA model.
Image Method Samp.Ratio PSNR (dB) Time (s) Rerror Brain 1 SBA 25.35% 34.63 5.32 0.067 512×512 Our model 25.35% 36.06 1.58 0.053  Brain 2 SBA 35.00% 40.04 8.47 0.050 512×512 Our model 35.00% 43.84 1.57 0.032  Brain 3 SBA 35.00% 36.34 1.44 0.029 200×200 Our model 35.00% 41.67 0.19 0.015Figure 3
2D brain image reconstruction with different sampling ratios. (a) The reference image; (b), (c), and (d) the recovered images with the sampling ratios of 7.6%, 15.0%, and 35.1%, respectively.
(a) (b) (c) (d)Figure4 shows that our model can also recover a faithful image with the sampling ratio of 35.0% of a real brain MR image with a higher level of complexity. The reference brain MR image is from the paper [30], and we denote it by Brain 2 in Tables 1 and 2.Figure 4
2D brain image reconstruction with different sampling ratios. (a) The reference image; (b), (c), and (d) the recovered images with the sampling ratios of 7.6%, 18.5%, and 35.0%, respectively.
(a) (b) (c) (d)In Figure5 we present the comparison results of our method and the SBA method in reconstructing three different real brain MR images, in which the image (i) can be found in [31] and denoted by Brain 3 in this paper. Column 1 shows the sample data, where image (a) is with the sampling ratio of 25.35% and images (e) and (i) are both with sampling ratios of 35.01%. Column 2 and Column 3 show the reconstruction results of the proposed model and the SBA model, respectively. The last column shows the reference images. We can see from Figure 5 that some fine structures of the reference image are missing with the SBA model, while our model can reconstruct clearer structures, which can be obviously observed by comparing (f) and (g).Figure 5
The comparison results of our method and the SBA method for three different brain MR images. Column 1: The sample data. Column 2: The reconstruction results of the proposed model. Column 3: The reconstruction results of the SBA model. Column 4: The reference brain MR images.In Figure6 we show the comparison results of our method and the SBA method with different sampling ratios in terms of the PSNR value and the relative error, where the results of the phantom image, the image Brain 1 and the image Brain 2 are shown in Column 1, Column 2, and Column 3, respectively. The first row of Figure 6 shows the PSNR values while the second row shows the relative errors. From Figure 6 we can see that PSNR values of the recovered images with our model are higher than those with the SBA model, and the relative errors are much smaller than those of the SBA model. This indicates that our model performs better in recovering MR images from undersampled images than the SBA model.Figure 6
The comparison results of our method and the SBA method with different sampling ratios. Row 1: The comparison of PSNR (dB) with the proposed model and the SBA model. Row 2: The comparison of relative errors with the proposed model and the SBA model.We evaluate the improvement of the proposed model by comparing the PSNR value and the computational time with the SBA model for a number of 100 brain MR images. The comparison results are shown in Figure7. Note that the sources of these brain images are the same as Brain 3 and the sizes of all images are 200×200. From Figure 7(a), we can observe that the PSNR values of reconstructed images with our model are all greater than those with the SBA model, which indicates the higher reconstruction accuracy of our model. Besides, from Figure 7(b), it can be seen the average computational time with our model is below 0.2s, much less than that with the SBA model, which is above 1.4s. We can conclude that our model is much more efficient and accurate than the SBA model.Figure 7
Comparison of reconstruction results between the proposed model and the SBA model with application to 100 brain MR images. (a) The PSNR (dB); (b) the computational time (s).
(a) (b)


## 4.3. Quantitative Analysis

In this section, to clearly see the numerical comparison of the reconstruction accuracy and efficiency with different sampling ratios for different images, we record the size, PSNR value, computational time, and relative error for the three images in Figures2, 3, and 4 in Table 1. From Table 1 we can observe that for the same image, with a higher sampling ratio, the PSNR value is larger and the Rerror value is smaller, which means that the reconstruction result is more accurate, while the CPU time is relatively less and thus is more efficient. This an intuitive and understandable result.We compare the numerical comparison results of our method and the SBA method for the three real brain MR images in Table2. Column 4 shows that the PSNR values of our model are 3.52 dB higher than those of the SBA model on average, which validates the improvement of reconstruction accuracy. Column 5 shows the computation time of the SBA model and the proposed model, from which we can obtain that the computation time of our model is much less than that of the SBA model. Besides, the Rerror values shown in Column 6 also verify higher reconstruction precision of our model than the other one. Table 2 demonstrates that our model can obtain reconstruction results with higher accuracy but less time. Therefore, our model is more accurate and efficient compared with the SBA method.


## 5. Conclusion

We present a fast compressed sensing model for accurate and efficient MR image reconstruction in this paper. By introducing a median filter into the proposed minimization problem, our model can reconstruct images and eliminate noise at the same time. Besides, because the proposed minimization problem is actually aL1 minimization problem, the split Bregman method can be applied to give a fast algorithm. We use lots of MR images to test the performance of our method and present reconstruction results. To quantitatively evaluate our model, we also give the comparison results between our method and the SBA method in terms of reconstruction accuracy and efficiency. Several indexes such as the PSNR values, the Rerror values, and the computation time are used to validate the performance of our model. It has been shown from experimental results that our model can more efficiently achieve higher precision of reconstruction results with the same sampling ratio. In summary, our model can reconstruct more accurate images with less computation time than the SBA model.



---

*Source: 8316194-2018-09-16.xml*