Mathematical Problems in Engineering

Volume 2018, Article ID 8316194, 9 pages

https://doi.org/10.1155/2018/8316194

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

^{1}School of Science, Harbin Institute of Technology, Shenzhen, China^{2}Department of Mathematics, Harbin Institute of Technology, Harbin, China

Correspondence should be addressed to Yunyun Yang; nc.ude.tih@nuynuygnay

Received 11 April 2018; Accepted 1 August 2018; Published 16 September 2018

Academic Editor: J.-C. Cortés

Copyright © 2018 Yunyun Yang 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

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.