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.
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 the 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.
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 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 is the measurement matrix, usually , is the measurement vector, and 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 -minimization problem:where the definition of is the number of nonzero elements in vector . 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 [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 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 aswhere represents the observed partial k-space data, represents the measurement matrix, and is a positive parameter. is the Fourier transform and 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].
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: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:where and are two positive parameters taking small values. is a two-dimensional median filter. Given an original image , the output image is defined aswhere is a template whose size is usually or .
The proposed minimization problem (8) is then converted into an unconstrained problem as follows:where 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 when the algorithm converges, where is the last iteration result. In our model, is replaced with and another TV norm is added into the problem (5). 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 will force to approach in each iteration, which means we can also achieve the effect of eliminating noise in the image . Therefore, by jointly minimizing the TV term , the data fitting term , and the median filter term , we can obtain a much clearer reconstructed image.
The median filter term 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 minimization problem; thus we can apply the split Bregman method to solve the problem (10). We first introduce two auxiliary variables and and obtain the following minimization problem:where and are two positive parameters. The two quadratic functions and are used to enforce weakly the constraints and . In order to establish strong constraints to the above two equations, we need to introduce additional relaxation variables and andthen the problem becomes the following sequence problems:where we update the Bregman variables and as follows:
We apply the alternating minimization scheme to solve the problem (12) with respect to , and as follows:
For the minimization problem (14), we can obtain the following equation with respect to :whereSince the Fourier transform is orthogonal, we have , let , 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
4. 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.
Algorithm 1.
Step 1. Input and , initialize , , , and , and give values for , , , and .
while do
Step 2. Update by (19).
Step 3. Update by (20).
Step 4. Update by (21).
Step 5. Update by (13).
end while
Step 6. Output the reconstruct result .
In Algorithm 1, the parameter is the tolerance criteria. The parameters , , , and 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 asFrom 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 aswhere is the reconstruction image and 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 aswhere 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.
(a)
(b)
(c)
(d)
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.
(a)
(b)
(c)
(d)
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.
(a)
(b)
(c)
(d)
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 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.
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 200200. 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.
(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 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.
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 a 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
The reference image (a) in Figure 2 is a synthetic image, which is available from the corresponding author upon request. 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]. 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]. 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]. 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.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by Shenzhen Fundamental Research Plan (no. JCYJ20160505175141489).