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

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 reconstructingMR 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.


Introduction
Compressed sensing [1][2][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][5][6][7], wireless sensor networks [8][9][10], speech coding [11,12], image memorability prediction [13], and others [14][15][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 timeconsuming.The split Bregman technique [17] is an efficient method which can be successfully used in image denoising, segmentation and compressed sensing [18][19][20][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 the  1 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 Section 2 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.

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  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: where  ∈ R × is the measurement matrix, usually  ≪ ,  ∈ R  is the measurement vector, and  ∈ R  is a sparse signal.Note that the measurement matrix  is usually chosen as the random matrices [22,23].
The reconstructed  can be realized by solving the following  0 -minimization problem: where the definition of ‖‖ 0 is the number of nonzero elements in vector .However, Natarajan [24] proved that the problem (2) is a nondeterministic polynomial hard (NPhard) problem in general.As an alternative, the convex relaxation of (2) has provided an approximation solution to  [25,26], which leads to the following convex optimization problem: Since the observed signal  usually contains noise, thereby the Lagrangian version of the problem (3) would be more suitable: where  is expected to the variance of noise.
The above  1 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 kspace 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.

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 where   represents the observed partial k-space data,  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: 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].

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

The Proposed Minimization Problem.
In this section, we present the following minimization problem: where  and  are two positive parameters taking small values. is a two-dimensional median filter.Given an original image V(, ), the output image (, ) is defined as where  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: 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 − V  ‖ → 0 when the algorithm converges, where V  is the last iteration result.In our model, V  is replaced with V  and another TV norm ‖V  − V‖ 1 is added into the problem (5).V  is the result by applying the median filtering  to the ℎ iteration reconstruction image, which can reduce the additional noise that comes from the reconstruction process, and ‖V  − V‖ 1 will force V to approach V  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 ‖V‖ 1 , the data fitting term ‖FV −   ‖ 2 2 , and the median filter term ‖V  − V‖ 1 , we can obtain a much clearer reconstructed image.
The median filter term ‖V  − V‖ 1 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.

Fast Minimization with the
We apply the alternating minimization scheme to solve the problem (12) with respect to V,  →  and  →  as follows: .
For the minimization problem ( 14), we can obtain the following equation with respect to V: where Since the Fourier transform is orthogonal, we have F T = F −1 , let  = ( T  − F T F + ), and then get For the other minimization problems ( 15) and ( 16), we can use the shrinkage operator to explicitly give the updating schemes for the variables  →  and  →  as follows: where the vector shrinkage operator is defined as

Numerical Algorithm and Experimental Results
4.1.Numerical Algorithm.Based on the above minimization schemes presented in Section 3.1, we give the following fast algorithm.

end while
Step 6.Output the reconstruct result V +1 .
In Algorithm 1, 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.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 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 where V  is the reconstruction image and V is the reference image.
In Figure 1 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 where sum(R) is the total number of the selected sampling pixels, while nembel(R) is the total number of pixels in R.
In Figure 2 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.
In Figure 3 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.
Figure 4 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.
In Figure 5 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).
In Figure 6 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     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.
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 Figure 7.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.

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 Figures 2, 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 Table 2. 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.

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 a  1 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.

Data Availability
(1) The reference image (a) in Figure 2 is a synthetic image, which is available from the corresponding author upon request.(2) The test image (a) in Figure 3 was used to support this study and is available at DOI: 10.1109/09/CVPR.2008.4587391.These prior studies (and datasets) are cited at relevant places within the text as [4].(3) The test image (a) in Figure 4 was used to support this study and is available at DOI: 10.1109/TMI.2010.2073717.These prior studies (and datasets) are cited at relevant places within the text as [30].(4) The test brain image (i) in Figure 5 and the 100 test brain images in Figure 7 were used to support this study and are available at http://www.engr.uconn.edu/∼cmli/.The images can be obtained by downloading the code attached to [31].(5) All other data, in Figures 1-7 and Tables 1-2, used to support the findings of this study, including the sample data, the reconstructed images, the error, the PSNR value, and the CPU computational time, are obtained by the authors implementing the SBA method and the proposed model for the test images.Requests for data will be considered by the corresponding author after publication of this article.
and give values for , , , and .

Figure 1 :
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.

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.
1 and the image Brain 2 are shown in Column 1, Column 2, and Column 3,

Figure 5 :
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.

Figure 6 :
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.

Figure 7 :
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).

Table 1 :
Numerical results for the synthetic phantom image and real MR images with different sampling ratios.

Table 2 :
Numerical comparison for different MR images between our model and the SBA model.