Medical image restoration method via multiple nonlocal prior constraints

Abstract

Medical image restoration is a fundamental issue in the area of medical signal processing, which aims recove high quality medical image from its degradation observation. Recently, the methods with nonlocal self-similarity prior have led to a great improvement on many medical image restoration tasks. Nevertheless, the nonlocal technique is generally embedded with only one kind of constraint, such as sparsity or low-rank in the conventional model, which limits their abilities and show good performance on certain prior. To address this problem, in this paper, we present a novel medical image restoration method with multiple nonlocal-based prior regularizations. The surfacelet transformation is introduced to construct a cubic sparsity constraint to a group of nonlocal similar patches. Likewise, due to the self-similarity existed in the medical image, two extra kinds of nonlocal-based priors, nonlocal total variation and nonlocal weighted low-rank, are also exploited to constrain the local smoothness and nonlocal relationship jointly. In this way, each of the designed priors can well recover a group of patches with similar structure. And then, the designed priors are combined into a unified proposed optimization framework, which will obtain the advantages from all of them simultaneously. Finally, to solve the objective function in the proposed framework, we develop an iterative numerical scenario based on alternating direction multipliers method. The extensive experiments on test medical images demonstrate that our proposed model outperforms the comparison methods on both of visual quality and objective evaluation results.

Keywords

Medical image restoration nonlocal prior low-rank sparsity surfacelet transformation

1 Introduction

MRI is an important method of medical diagnosis, which is accurate and has little damage to the human body. However, MRI also has the disadvantage of excessive scanning time, resulting in motion artifacts due to the movement of the patient. And it is not conducive to obtaining correct examination results. Besides, another challenge in MRI processing is the noises produced by the natural attribute of the imaging system, the motion displacement of the patent, the existence of some metal foreign body, the influence of the external environment and etc.

Researchers have proposed a flurry of methods to solve these problems, including accelerating the sampling procedure, and having the origin image restored. The reconstruction and denoising of MRI image can be treated as part of image restoration task and be both modeled by the following degenerate process. $y = Hx + n$ (1) where x ∈ i^N is the original image, y ∈ i^M is our measurement, n ∈ i^M×1 is some additive noise, H ∈ i^M×N is a known measurement/sensing matrix and its form depends on different specific restoration task. Considering the base case, if H is matrix of full column and n is a zero vector. It means x is simple linear transformed to y without losing any information.

Then the task of recovering x is completed by solving the following equation $min_{x} ∥ y - Hx ∥$ . The solution becomes x = (H^TH) ^-1H^Ty.

And the solution is more complicated in other situations. If H is a undersampling matrix in some transform domain, which always happen when the signal is sampled below the Nyquist rate, it’s a highly underdetermined problem to reconstruct x from undersampled data y because that there are more variables than the equation. There exist lots of possible solutions. However, if x is sparse, x can be reconstruct in a high possibility according to the CS theory [1, 2] with a measurement matrix H that is highly incoherent with Φ. x is sparse if it can be written as x = Φc, where Φ = {Φ₁, Φ₂, ⋯ , Φ_N} is so called sparsifying transform domain basis of i^N and c is a vector containing much more zero entries than nonzero entries, written as ∥c ∥ _{l
₀} = k = N.

According the CS theory, the estimation of the restoration task is reconstructed by linear programing $min {{∥ x ∥}_{l_{1}} : y = Hx}$ (2)

In order to gain better robustness to noise, different kinds of regularization prior (i.e., sparse prior) about the image are added as penalty term R (x). $min_{x} {\frac{1}{2} ∥ y = Hx ∥ + λ R (x)}$ (3)

The first term punishs those solution based on distance measure (i.e., mean square error ${∥ y - Hx ∥}_{2}^{2}$ ), to restrict the deviation from measurement. The second term can be one or combination of more penalty term x. And the parameter λ control a trade-off between two terms. Different penalty term is developed to apply different image prior. For instance, in image denoising task, the matrix H is an identity matrix, and n is an additive noise like Gaussian white noise, or burst noise and etc. And the actual probability distribution of n has is determined by the degenerate process. Usually the estimate is similar to Equation (2), but the second regularization term can be sparsity regularization term (i.e., ∥x ∥ _{l
_p}, ∥Φx ∥ _{l
_p}) Total variation (TV) norm [3], wavelet, Bregman distance [4] or low-rank. There are two popular choices for the discrete TV.

The first is isotropic and the second one is anisotropic l₁-based TV [5]. Total variation has been introduced by Rudin-Osher and Fatemi (ROF) in [3] as an edge-preserved local image transform. Because the image polluted by the noise have a larger total variation than noise-free image, it can indicate the smooth degree and is used widely in deburring and denoising.

Besides TV norm, low rank is another popular approach that capture the nature of the signal like images and videos. Hoverer, minimizing rank (x)is proved a NP hard problem, so we minimizing the convex envelope of rank (x), so-called nuclear norm instead.

A low-rank matrix refers to a matrix with a linear correlation between rows and columns or columns of a matrix, which has generally only a few large singular values. With the low-rank matrix completion, some elements of a low-rank matrix lost or damaged can be recovered according to the known elements. The training sample is divided into multiple sub-blocks, then the low-rank matrix recovery algorithm is used to learn from the sub-blocks and find the weight of the image used in block reconstruction. Lingala et al. [6] have also found a way to apply both low-rank constraints and sparse priors to the reconstruction of magnetic resonance images. And high acceleration and background separation are achieved by apply low-rank in a decomposition model [7, 8].

In addition to exploiting a prior information about original image, developing effective methods to exploit the self-similarity has been the topic of much research recently. Non-local methods assume that structural self-similarity in different patches of an image is ubiquity. Therefore, redundant information existing in images patch with similar texture and structures can restrain noise and provide detail for each other. And nonlocal regularizers have shown significant improvement on performance in image restoration tasks. Enlighten by successful application of nonlocal means, several non-local methods have also been proposed in image filtering [9], denoising [10, 11] and dehazing [12].

Due to the advantage of the nonlocal technique in image restoration, we present a novel approach that designs three nonlocal embedded constraint and puts them into a unified optimization framework, which can not only obtain better visual quality but also preserves more fine structures, to address the medical image restoration problem. The main contributions of our work are as follows.

We design a nonlocal sparsity constraint based on surfacelet transformation. To exploit the redundant information, we firstly group the self-similarity nonlocal patches into cubic matrix as a basic processing unit. And then, we constrain the sparsity of the unit by surfacelet transformation, which can better preserve the two dimensional structures of image patches.

We combine the designed cubic nonlocal sparsity constraint, the weighted nonlocal low-rank constraint and nonlocal total variation into a unified framework, which can exploit the advantages of the three regularization items simultaneously.

An iterative numerical scheme based on alternating direction multiplier method is developed to solve the objective function in the proposed framework efficiently and effectively.

At last, we implement extensive test experiments and compare our proposed with some existing excellent restoration algorithms on two representative restoration tasks, image compressive sensing and denoising. The results of the experiments show that our method outperforms others for various noise levels and compressive sampling rates.

2 Related work

Nowadays many methods have been proposed for image processing. Such as image denoising, image reconstruction [13–19 , 20] and basic theories of image processing [21 –23]. This method can be used in the people’s life, such in the medical area and so on.

The work in [24] proposed a way to deal with the image denoising, through this way can get a computationally simple denoising method by using the LRA. This method is different from the other SVD-based methods, the LRA in SVD domain avoids learning the basis for representing image patches, which usually is computationally expensive. The result of this method shows that it can effectively reduce noise in the image denoising. Besides this way of imaging denoising, in work [25] also proposed an efficient way, they proposed an implementation of distributed parallel optimization of K-SVD method on Spark, the result of the experiment shows that this method can not only get a good speed-up ratio, but also remains the details of the image. In work [26], they proposed a new method of image denoising, which can be used for completion and denoising of multilinear data.

Besides image denoising, in order to get better image, image reconstruction is also very important. Recent years, there are lots of new methods have been proposed for image reconstruction, in order to get a better result of the image reconstruction, there are a lot of papers which proposed new methods to obtain better results. As for the recovery of the image, Bayesian methods can be used, in work [27] they use the Bayesian method to recovery the images which are incomplete. The work in [14] W. Dong, G. Shi, X. Wu, and L. Zhang proposed use the AR models to construct adaptive sparsity regularizations for CS image recovery. In work [15], they proposed a new method to make the image restoration which is named SAIST. In this algorithm, they take a low-rank method which towards SSC and provide a conceptually simple interpretation from a bilateral variance estimation perspective. They proposed a new class of image recovery which is called SAIST. The work in [25], they proposed a new method of image restoration called SAIST. This algorithm generalizes the Bayes-Shrink from local to nonlocal models in order to solve the noise data. After the experiment, the result shows that this method can solve the image processing efficiently, especially when the noise levels is high and in the situation with large amount of missing data.

The basic theory of image reconstruction is also important, as for the theory in [28] proposed a novel iteratively reweighted TV regularization for CS reconstruction. This method is mainly used the posteriori estimation of the image gradients to compute the weight of spatially adaptive. Yue Lu and Minh N. Do combined the multi-scale pyramid and multi-directional filter bank (NDFB) and proposed the Surfacelet transform [29], which is suitable for the processing of 3D images. They proposed a new family of filter banks, named NDFB, which can achieve the structured construction. In work [30] presents a novel image classification framework by using the low-rank matrix and Laplacian group sparse. In order to improve the deficiency of the total variation in compressive sensing, they proposed a new method in work [13], which can solve a combination of total variation and a new regularization term under the framework of compressive sensing.

Total variation (TV) is a method which uses the fact for most of the images have a detail distribution of gradients, for that TV [31 –33] can be used as a model which can be modeled as a Laplace distribution. In work [31] proposed a new method of CS, for that by using an adaptive reweighted TV strategy to get a much better preserve image. They also get the weight of the low rank regularization by using the redundancy of non-local image patches. The result of the experiment shows that their method is much better than the others approaches. The work in [32] proposed a new compressed sensing method which combines both of the external and the internal information together to get the high-performance reconstruction of image reconstruction. They use the probabilistic atlas to control the level of the gradient regularization at each image location, within a weighted TV regularization prior. The result of the experiments show that this method is better than the other approaches, especially for different sampling rates and noise levels. The work in [33] proposed a method of saving the gradient histogram. They invented an efficient optimization method, which was developed and has achieved a satisfied result for kernel estimation. They combine the gradient histogram preservation prior with the previous method of image deburring, the new method shows that in this way can improve the simulations. And the result of the new method shows that it can achieve a high SSIM.

Recently, low rank approximation [34], low rank representation [35], sparse representation [36] have shown strong capabilities in signal approximation and subspace separation, which attract attention of the academic circles recently. In work [37], they made a very comprehensive theoretical analysis of low rank matrix recovery, and proposed a matrix decomposition into the sum of low rank matrices and sparse matrices. The representation-based theory has been widely used in artificial intelligence, image processing, pattern recognition [38, 39], computer vision and other fields.

3 Our proposed method

With the deep analysis about the previous related works, we can see that the nonlocal technique is great effective for image restoration due to the redundant local similar structures included in the image itself. And also, some ripe applications of nonlocal technique has been realized in many practical visual restoration tasks. Thanks for the good ability in preserving the textures and edges, the nonlocal technique based methods achieved excellent results and attracted more and more attentions recently. In light of the advantages of nonlocal technique, we proposed three kinds of nonlocal-based regularization constraints, which embedded the nonlocal techniques into the conventional low-rank, total variation and sparse models. And then, a comprehensive framework will be established with the mentioned three kinds of regularizations jointly to address the medical image restoration problems.

3.1 Weighted nonlocal low-rank constraint

Low rank prior attracted many attentions recently due to its good performance in area of signal processing. To our best knowledge, the low rank model, which is first proposed in robust PCA, can find the latent subspace structure of signals and extract the relationship among them. With the robust PCA, the following low rank model is implemented on different images to address background modeling in [40] as $\begin{matrix} min {∥ X ∥}_{*} + {∥ E ∥}_{2, 1} \\ s . t . Y = X + E \end{matrix}$ (4) where Y is a matrix, which includes the continuous frames of video as its column. X denotes the estimated recovery component, and E denotes the error matrix to obtain the outliers existed in these images. With the low rank minimization constraint on X, the correlation information can be exploited effectively. For these frames, the correlation component is the fixed background images, and the outlier components among them are just the moving foreground components.

Though the successful applications in background subtraction has been achieved with Equation (4), for single medical image, there is not enough related images that can be used for restoration problem. Hence, our proposed method adopted a patch-based prior model to restore the medical image from its compressive observation or degradation version, which can search some related patches to help us to establish the low rank constraint. In fact, similar to the natural image, there are a lot of self-similarity patches existed in the nonlocal areas of medical image itself. So, different from the traditional low rank model for video, we can impose the low rank constraint on the groups of similar patches in single image. Some studies have shown that the patches with similar structures lie in the same low-dimensional subspace and the matrix composed of these patches has low rank property.

Let x_i ∈ Rⁿ be a patch image with the size of $\sqrt{n} \times \sqrt{n}$ centered around the i-th pixel in image x, which is expressed as x_i = R_ix and R_i denotes a patch selection matrix. Then, for each pixel i, we define a patch group P_i ∈ R^n×N to include the similar patch $R_{i}^{k} x$ , where $R_{i}^{k}$ denotes selection matrix to extract the k-th most similar patch of x_i and N is the predefined total number of similar patches. Denote $x_{i}^{k}$ as $x_{i}^{k} = R_{i}^{k} x$ , we can obtain the similar patch group as $P_{i} = [x_{i}^{1}, x_{i}^{2}, \dots, x_{i}^{k}, \dots, x_{i}^{N}]$ . For the similar measurement, we adopt the following Euclidean distance as ${∥ x_{i}^{k} - x_{i} ∥}_{2} < d_{0}$ (5) where d₀ is a user-defined distance constant to constrain the upper boundary of correlation among patches.

With the matrix P_i, we constructed a nonlocal low rank prior with the weighted nuclear norm in [41] as $LR (x) = \sum_{i} {∥ P_{i} ∥}_{*, ω} = \sum_{i} \sum_{t} ω_{t} σ_{t} (P_{i})$ (6) where the ω = [ω₁, ω₂, ⋯ , ω_t] is the weighted coefficient vector, which satisfies ω_t ⩾ 0. ω_t denotes a nonnegative coefficient assigned to the t-th singular value of matrix P_i. With the Equation (6) the weight vector will enhance the representation capability of the original nuclear norm. Rational weights specified based on the prior knowledge and understanding of the problem will benefit the corresponding weighted nuclear norm minimization model for achieving a better estimation of the latent structures from the degradation data. In addition, each member of the matrix P_i is similar to each other in structures, and the low-rank based nuclear norm will help to recovery the latent related structures hidden in the nonlocal data group. Hence, take the nonlocal similar patches as the processing unit will benefit the designed prior in Equation (6) and obtain better visual recovery quality.

3.2 Nonlocal total variation

Introduced by Rudin-Osher and Fatemi (ROF) in [3], total variation is a noising removing algorithm that has superb performance comparing with least square methods based on L₂ norm. And it can provide sharp image without spurious noise. Moreover, the solving of constrained minimization problem is fast and simple. TV has been developed and combined with other methods [42 –44] since it is proposed by ROF. $R (u) = \int_{Ω} | \nabla u | d x d y$ (7)

Regularization term proposed by ROF is basically L₁ norm of the gradient.

where $| \nabla u | = \sqrt{u_{x}^{2} + u_{y}^{2}}$ for u = u (x, y) is the expected intensity function of a clean image and Ω ∈ i² is a bounded open set in i².

In order to acquire a better presentation of the image, non-local TV is used in our model. Comparing to NWTV, NLTV [45] does not weight the gradient by the magnitude of the gradient. Instead, the weight is determined not only by local structural around x but also the spatial similarity of neighborhood area of image.

Compared with TV, the non-local penalty term considers both the position information and gray level information of all pixels in the neighborhood of the pixel. Therefore, R_NLTV can describe the texture and detail adaptively according to patch-by-patch structure similarity, and the edge and texture is preserved better in reconstruction of image.

In view of the advantage of nonlocal total variation regularization, we introduced the NLTV as one of the regularization members of our proposed nonlocal-induced framework for medical image restoration. In the conventional total variation framework, the recovered image will generate much artifacts and serious stairs effect, which would loss many details of texture structures and edges. Nevertheless, by embedding the nonlocal technique, it is believed that more details and local textures would be preserved to improve the visual quality.

3.3 Nonlocal sparsity constraint

The similar image patches are stacked into a three- dimensional array in our model, and the image stack can be seen as a 3D signal. Because of the similarity of structural and textual among these patches, and correlation among the pixels in a patch, the 3D signal has high redundancy. And this fact suggests us that with proper transform, an efficient representation for the signal is possible. And if the redundancy is exploited in denoising, the reconstruction algorithm can achieve better robust to the noise.

Surfacelet Transform has a strong ability for signal processing. It has strong multi-dimensional signal processing capability and implements multi-directional multi-scale decomposition of multi-dimensional signals. It has been used in video denoising application [46, 47] for its characteristic of low redundancy and high computation efficiency. The video is a special 3D signal with one temporal dimension and two spatial dimensions. And the 3D signal formed by similar images patches have also three dimensions.

In light of the advance of surfacelet, we collect a cubic matrix with the similar nonlocal patches extracted from the original medical image, and then enforce the sparsity constraint on it with surfacelet transformation. The mentioned sparsity constraint S (x) is designed as follows: $S (x) = \sum_{i} {∥ ST (M (x_{i})) ∥}_{1}$ (8) where M (x_i) is the collected cubic matrix, which includes the nonlocal similar patches of the i-th patch. ST (•) denotes the surfacelet transformation that can compute the transformation coefficiets of the processed unit M (x_i), and ‖‖₁ denotes the l₁ norm, which can constrain the sparsity of coefficient. By the Equation (8), the patches with similar structures are processed as a whole unit, and it is believed that the designed regularization will preserve more details in the recovered image. Also, different from the conventional sparsity-based nonlocal technique, our proposed constraint processes the unit in a cubic way, which will help to preserve more two dimensional structure existed in image patch effectively.

3.4 The proposed comprehensive framework

In this section, we give our proposed objective function by combining the above three nonlocal-based constraints as follows. $\begin{matrix} min_{x} \frac{1}{2} ∥ y - Hx ∥ + λ_{1} R_{NLTV} (x) \\ + λ_{2} LR (x) + λ_{3} S (x) \end{matrix}$ (9) where the last three terms are the constraint detailed in the above subsections. λ₁, λ₂, λ₃ are the positive scalar to balance them. With the proposed objective function in Equation (9), we give a nonlocal-based framework for medical image restoration problem. The comprehensive framework combines three nonlocal-based regularization prior constraints to benefit from all of them simultaneously. The second term constrain the nonlocal smoothness of the estimated image. The third term aims to recover the related structures hidden in the similar patches by exploring their low rank component. To preserve the structure from the view of two dimension, the last term arranges a 3D array with similar patches, and enforce the sparsity of ST coefficient due to the redundant structural information among similar patches. It is noted that the three constraints will guide the optimization to find a solution to satisfy all of the requirements, which can earn different recovered property and preserve more fine structures to improve visual quality greatly.

4 The solution scheme for the proposed objective function

In this section, we will present a efficient numerical solution to solve the proposed objective function in Equation (9) Due to the three nonlocal-based constraints, solving the problem in Equation (9) is a challenging task. To obtain a high quality recovered image, we adopt the iterative numerical scheme based on alternating direction multiplier method (ADMM) to get an approximate convergence solution.

To decouple the variables, we introduce three auxiliary variables u, v, q and then the optimization in Equation (9) can be convert into the following constraint problem: $\begin{matrix} min_{x} \frac{1}{2} ∥ y - Hx ∥ + λ_{1} R_{NLTV} (x) \\ + λ_{2} LR (x) + λ_{3} S (x) \\ s . t . u = x, v = x, q = x \end{matrix}$ (10)

And then, the augmented Lagrange function of constrained Equation (10) can be obtained as follows. $\begin{matrix} min_{x} \frac{1}{2} ∥ y - Hx ∥ + λ_{1} R_{NLTV} (x) \\ + λ_{2} LR (x) + λ_{3} S (x) + \frac{μ_{1}}{2} {∥ u - x ∥}_{2}^{2} \\ + 〈 u - x, θ_{1} 〉 + \frac{μ_{2}}{2} {∥ v - x ∥}_{2}^{2} \\ + 〈 v - x, θ_{2} 〉 + \frac{μ_{3}}{2} {∥ q - x ∥}_{2}^{2} \\ + 〈 q - x, θ_{3} 〉 \end{matrix}$ (11) where μ₁, μ₂ and μ₃ are three positive regularization parameters to balance the terms. θ₁, θ₂ and θ₃ are the Lagrange multipliers corresponding to each auxiliary variable. To solve it easily, a more compact equivalent form is shown as $\begin{matrix} min_{x, u, v, q} \frac{1}{2} ∥ y - Hx ∥ + λ_{1} R_{NLTV} (u) \\ + λ_{2} LR (v) + λ_{3} S (q) + {∥ u - x + \frac{θ_{1}}{μ_{1}} ∥}_{F}^{2} \\ + {∥ v - x + \frac{θ_{2}}{μ_{2}} ∥}_{F}^{2} + {∥ q - x + \frac{θ_{3}}{μ_{3}} ∥}_{F}^{2} \end{matrix}$ (12)

Next, with the ADMM, we can solve the optimization in Equation (12) iteratively. At the k-th iteration, by fixing other variables, we can firstly obtain an x subproblem about as $\begin{matrix} min_{x} \frac{1}{2} {∥ y - Hx ∥}_{2}^{2} + {∥ u^{k} - x + \frac{θ_{1}^{k}}{μ_{1}} ∥}_{F}^{2} \\ + {∥ v - x + \frac{θ_{2}^{k}}{μ_{2}} ∥}_{F}^{2} + {∥ q - x + \frac{θ_{3}^{k}}{μ_{3}} ∥}_{F}^{2} \end{matrix}$ (13)

Equation (13) is a typical convex quadratic cost function which can be solved easily, and its closed-form solution can be achieved by forcing its derivate to be zero. $\begin{matrix} x^{k + 1} = 2 * {(H^{T} H + 6 I)}^{- 1} \\ * (\begin{matrix} (u^{k} - \frac{θ_{1}^{k}}{μ_{1}}) + (v^{k} - \frac{θ_{2}^{k}}{μ_{2}}) + \\ (q^{k} - \frac{θ_{3}^{k}}{μ_{3}}) + H^{T} y \end{matrix}) \end{matrix}$ (14)

Removing the terms independent of u, we have $min_{u} {∥ u - x^{k + 1} + \frac{θ_{1}}{μ_{1}} ∥}_{F}^{2} + λ_{1} R_{NLTV} (u)$ (15)

The minimization in (15) is the so-called nonlocal total variation recovery problem. Its solution can be computed with the method presented in [45].

Similarly, by removing terms independent of v, we can obtain $min_{v} {∥ v - x^{k + 1} + \frac{θ_{2}}{μ_{2}} ∥}_{F}^{2} + λ_{2} LR (v)$ (16)

Incorporating Equation (10) into Equation (16), the minimization can be converted into the following form $min_{v} {∥ v - x^{k + 1} + \frac{θ_{2}}{μ_{2}} ∥}_{F}^{2} + λ_{2} \sum_{i} {∥ P_{i}^{v} ∥}_{*, ω}$ (17)

Assume that $Q_{i}^{k + 1}$ is a grouping matrix from $x^{k + 1} - \frac{θ_{2}}{μ_{2}}$ , and $Q_{i}^{k + 1} = U_{Q_{i}} Σ_{Q_{i}} V_{Q_{i}}^{T}$ . For the above optimization in Equation (17), it can be solved by $V_{i}^{k + 1} = U_{Q_{i}} Shrinkage (Σ_{Q_{i}}) V_{Q_{i}}^{T}$ (18) where shrinkage (∑_{Q
_i}) is the soft-shrinkage operator and its diagonal element is computed as $max ({(Σ_{Q_{i}})}_{jj} - ω_{j}^{i}, 0)$ . Then, we can divide the grouping patch with the inverse operator to aggregate $v_{i}^{k + 1}$ .

According to the reweighted scenario presented in [48], we define the j-th element of weighted vector ω as follows: $ω_{j}^{i} = \frac{m}{σ_{j} (P_{i}) + ɛ}$ (19) where m is a constant and ɛ is a small positive scalar to avoid zero value. σ_j (P_i) is the j-th singular value of P_i.

Next, we obtain the subproblem for q by fixing others as $min_{q} {∥ q - x^{k + 1} + \frac{θ_{3}}{μ_{3}} ∥}_{F}^{2} + λ_{3} S (q)$ (20)

With the definition in Equation (20), q is solved by $q = {ST}^{- 1} ({Shrinkage}_{λ_{3}} (ST (M (q_{i}))))$ (21)

With the above description, our comprehensive proposed algorithm is summarized in detail in Algorithm 1.

Algorithm 1
Input: The observed image y, degradation kernel H, maximum
iterative number K;
Output: The restoration image x;
While not convergence do
1. Compute the x^k+1 with Equation (14);
2. Compute the u^k+1 with Equation (15) by the method in [45];
3. For each i, compute the $V_{i}^{k + 1}$ with Equation (18);
4. Compute the v^k+1 by aggregating ${V_{i}^{k + 1}}_{i = 1, 2 . . .}$
5. Update the parameters:
$θ_{1}^{k + 1} = θ_{1}^{k} - μ_{1} (x^{k + 1} - u^{k + 1})$
$θ_{2}^{k + 1} = θ_{2}^{k} - μ_{2} (x^{k + 1} - v^{k + 1})$
$θ_{3}^{k + 1} = θ_{3}^{k} - μ_{3} (x^{k + 1} - q^{k + 1})$
6. k = k+1;
End

5 Experiments and analysis

In this section, we evaluate the performance of the proposed method and compare it with several benchmark methods. The performances of our experiments are evaluated on four gray medical images. The test images are shown in Fig. 1. We perform experiments on two kinds of restoration tasks that are compressive imaging and denoising. Due to the limited space, we only show parts of the visual restoration results. To evaluate the performance objectively, we take the PSNR as the objective evaluation indicator.

Fig. 1

Medical images for testing.

5.1 Experiments on compressive imaging

When testing the performance on medical images compressive reconstruction, the degradation kernel H is defined as an under-sampling matrix, which only selects parts of the original data for generating the observation image. Our proposed method is compared to some recent excellent algorithms for compressive reconstruction. The comparison methods include total variation-based nonlocal compressive reconstruction method (TVNCR) in [49], the collaborative sparsity-based compressive reconstruction method (CSCR) in [50], block compressive reconstruction with directional transforms (BCRDT) [51], compressive reconstruction with multi-hypothesis (CRMH) [52] and nonlocal low-rank based compressive reconstruction (NLRCR) proposed by [53]. The parameters for our proposed method are set empirically as: λ₁ = 2 . 0, λ₂ = λ₃ = 1.0, μ₁ = μ₂ = 0.1, the maximum iterative step K = 80. The patch size in our experiment is set to be 8 ×8, and the number of nonlocal similar patch used for designed regularization is set to be 20. The parameters for other comparison methods are set to be same with their papers. The experiments are implemented six times on each test image and the average PSNR result versus undersampling ratio is reported finally in Fig. 4. Due to the limited space, we only show part of the visual reconstructed results in Figs. 2–4 when the undersampling ratio is 60%.

Fig. 2

The visual reconstruction results on image 1.

Fig. 3

The visual reconstruction results on image 2.

Fig. 4

The PSNR curve of reconstruction result by comparison methods.

From the visual results, it can be seen that our proposed method shows better performance on edge and details recovery, which is significant quality evaluation of restoration tasks. Additionally, our method also obtains the highest objective indicator among the comparison methods, and shows obvious advantage when under-sampling ratio is low. In a word, both of the visual and objective results verify the competitive performance and effectiveness of our proposed method.

5.2 Experiments on denoising

In this subsection, we test the denoising performance on medical images. In denoising experiment, the degradation kernel H is defined as an identity matrix, and extra noise is introduced into the observation image. In our experiments, the Gaussian white noise is used and the standard deviation is varied from 5 to 35 with 10 as interval. The compared algorithms include total variation-based denoising method called TVAL in [54], Image denoising via simultaneous sparse coding (IDSSC) in [55], nonlocal joint regularization-based image denoising method (IDNJR) [56], image denoising with gradient histogram preservation (IDGHP) [57 –60], and weighted nuclear norm-based image denoising (IDWNN) method in [21]. The parameters for our proposed method are set same as the compressive reconstruction experiments. Similarly, we only show the part of denoising visual results when the noise deviation is 20 in Figs. 5, 6 respectively. The curves of average PSNR and SSIM results versus noise standard deviation are plotted in Figs. 7 and 8 to evaluate the denoising performance objectively.

Fig. 5

The denoising visual results on image 1.

Fig. 6

The denoising visual results on image 3.

Fig. 7

The PSNR results of test images.

Fig. 8

The SSIM results of test images.

From the visual results, all the comparison methods remove the most of noise existed in medical image. However, our method shows advantage on fine structures preservation and improve the visual quality greatly. The highest PSNR and SSIM curves also demonstrate that our proposed denoising method outperforms other existing excellent methods.

6 Conclusion

In this paper, a novel nonlocal-based multiple constraints framework is presented to address the medical images restoration problem. The proposed method exploits three kinds of patch-based nonlocal priors, which are nonlocal weighted low-rank, nonlocal total variation and nonlocal sparsity under surfacelet transformation, as constraints and incorporated them into a unified framework to gain the advantage from these priors simultaneously. To further improve the performance on details recovery, we introduce the surfacelet transformation to construct a cubic sparsity prior, which can better preserve the two-dimensional structures of nonlocal patches. Different from the conventional approaches, which focus on either nonlocal self-similarity or local smoothness preservation, our proposed method jointly constraint the two aspects in the same framework. Furthermore, an iterative optimization technique based on ADMM is proposed to solve the framework effectively to obtain an approximate solution.

The performance of our method is evaluated on two kinds of representative medical image restoration tasks, medical image compressive reconstruction and denoising. Various visual recovered results show the advantage of our proposed method on fine structures preservation. And also, the PSNR curves verifies the superiority compared to current conventional algorithms.

References

Donoho

D.L.

, Compressed sensing[J], IEEE Transactions on Information Theory 52(4) (2006), 1289–1306.

Lustig

, Donoho

and Pauly

J.M.

, Sparse MRI:, The application of compressed sensing for rapid MR imaging[J], Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 58(6) (2007), 1182–1195.

Rudin

L.I.

, Osher

and Fatemi

, Nonlinear total variation based noise removal algorithms, Phys D, Nonlinear Phenomena 60(1–4) (1992), 259–268.

Bregman

L.M.

, The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming, Comput Math Math Phys 7(3) (1967), 200–217.

Beck

and Teboulle

, Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[J], IEEE Transactions on Image Processing 18(11) (2009), 2419–2434.

Otazo

, Candès

and Sodickson

D.K.

, Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components[J], Magnetic Resonance in Medicine 73(3) (2015), 1125–1136.

Jin

K.H.

, Lee

and Ye

J.C.

, A general framework for compressed sensing and parallel MRI using annihilating filter based low-rank Hankel matrix[J], IEEE Transactions on Computational Imaging 2(4) (2016), 480–495.

Haldar

J.P.

, Low-Rank Modeling of Local $ k $-Space Neighborhoods (LORAKS) for Constrained MRI[J], IEEE Transactions on Medical Imaging 33(3) (2014), 668–681.

Zhang

, Cao

and Wang

, A new image filtering method: Nonlocal image guided averaging[C], Acoustics, Speech and Signal Processing (ICASSP), 2014 IEEE International Conference on IEEE (2014), 2460–2464.

10.

Jia

, Feng

and Wang

, Rank constrained nuclear norm minimization with application to image denoising, Signal Processing 129 (2016).

11.

Chen

, Tang

, Zhang

and Lei

, Image decomposition and denoising based on shearlet and nonlocal data fidelity term, Signal Image & Video Processing (7) (2018), 1–8.

12.

Berman

and Avidan

, Non-local image dehazing[C], Proceedings of the IEEE conference on computer vision and pattern recognition (2016), 1674–1682.

13.

Canh

T.N.

, Total variation reconstruction for kronecker compressive sensing with a new regularization, in: Picture Coding Symposium (2013), pp. 261–264.

14.

Dong

, Shi

, Wu

and Zhang

, A learning-based method for compressive image recovery, J Vis Commun Image Represent 24(7) (2013), 1055–1063.

15.

Dong

, Shi

and Li

, Nonlocal image restoration with bilateral variance estimation: A low-rank approach, IEEE Trans Image Process 22(2) (2013), 700–711.

16.

, Chen

, Lin

and Sun

, Hyperspectral Image Denoising with Composite Regularization Models.[J], Journal of Sensors 2016.

17.

and Sun

, Image completion approaches using the statistics of similar patches, IEEE Trans Pattern Anal Mach Intell 36(12) (2014), 2423–2435.

18.

Zhang

, Desrosiers

, Qu

, Guo

and Zhang

, Medical image super-resolution with non-local embedding sparse representation and improved IBP, in Proc IEEE Int Conf Acoust, Speech Signal Process. (ICASSP) (2016), pp. 888–892.

19.

Dong

, Zhang

, Shi

and Li

, Nonlocally centralized sparse representation for image restoration, IEEE.

20.

and Qi

, A douglas-rachford splitting approach to compressed sensing image recovery using low-rank regularization, IEEE Transactions on Image Processing 24(11) (2015), 4240–4249.

21.

, Xie

, Meng

, Zuo

, Feng

and Zhang

, Weighted nuclear norm minimization and its applications to low level vision, Int J Comput Vis 121(2) (2016), pp. 183–208.

22.

Zhang

, Matuszewski

B.J.

, Shark

L.K.

and More

C.J.

, Medical image segmentation using new hybrid level set method, International Conference Biomedical Visualization 2008.

23.

Leo Alexander

Dr. and Monica

, Statistical and Machine. Learning Techniques for Prediction of Customer Churn in Telecom, International Journal of Innovative Research in Computer and Communication Engineering (IJIRCCE), V5, I5, 2017.

24.

Paul

, Wu

, Yang

J.F.

and Jeong

, Gradient-based edge detection for motion estimation in H. 264/AVC, IET Image Processing 5(4) (2011), 323–327.

25.

Guo

, Zhang

and Liu

, An efficient SVD-based method for image denoising, IEEE Trans Circuits Syst Video Technol 26(5) (2016), 868–880.

26.

Tan

, Wei

, Wu

Z.B.

, Chen

and Gu

, Parallel optimization of ksvd algorithm for image denoising based on spark, in: IEEE International Conference on Signal Processing (2017), pp. 820–825.

27.

Zhang

, Ely

, Aeron

, Hao

and Kilmer

, Novel methods for multilinear data completion and de-noising based on tensor SVD, in Proc IEEE Conf ComputVisPattern Recognit (2014), pp. 3842–3849.

28.

Zhou

, et al., Nonparametric Bayesian dictionary learning for analysis of noisy and incomplete images, IEEE Trans Image Process 21(1) (2012), 130–144.

29.

Dong

, Yang

and Shi

, Compressive sensing via reweighted TV and nonlocal sparsity regularisation, IET Electron Lett 49(3) (2013), 184–186.

30.

and Do

M.N.

, Multidimensional directional filter banks and surfacelets. [J], IEEE Trans Image Processing 16(4) (2007), 918–931.

31.

Zhang

and Ma

, Low-rank decomposition and laplacian group sparse coding for image classification, Neurocomputing 135 (2014), 339–347.

32.

Zhang

, Desrosiers

and Zhang

, Effective compressive sensing via reweighted total variation and weighted nuclear norm regularization, Acoustics, Speech and Signal Processing (ICASSP), 2017 IEEE International Conference on IEEE (2017), pp. 1802–1806.

33.

Zhang

, Desrosiers

and Zhang

, Atlas-based reconstruction of high-performance brain MR data, Pattern Recognit 76 (2018), 549–559.

34.

Zhang

, Desrosiers

, Zhang

and Cheriet

, Effective document image de- blurring via gradient histogram preservation, Image Processing (ICIP), 2015 IEEE International Conference on. IEEE, (2015), pp. 779–783.

35.

Markovsky

, Low rank approximation: Algorithms, implementation, applications[M], Springer Science and Business Media, 2011.

36.

Liu

, Lin

and Yu

, Robust subspace segmentation by low-rank representation[C], Proceedings of the 27th Conference on Machine Learning (ICML-10) (2010), 663–670.

37.

Huang

, Zhang

and Metaxas

, Efficient mr image reconstruction for compressed mr imaging, Medical Image Analysis 15(5) (2011), 670–679.

38.

Candes

E.J.

, Li

, Ma

, et al., Robust principal component analysis? [J], Journal of the ACM(JACM) 58(3) (2011), 11.

39.

, Wu

, Chen

, et al., Collaborative Self-Regression Method With NonlinearFeature Based on Multi-Task Learning for Image Classification.[J] IEEE Access 6, pp. 43513–43525.

40.

, Chen

, Wu

, et al., Self-supervised sparse coding scheme for imageclassification based on low rank representation.[J] PLOS One 13(6) (2018), 1–15.

41.

, Wang

and Liu

, Low-Rank and Sparse Structure Pursuit via Alternating Minimization, in Proc of the 19th International Conference on Artificial Intelligence and Statistics, 2016.

42.

, Xie

, Meng

, Zuo

, Feng

and Zhang

, Weighted nuclear norm minimization and its applications to low level vision, International Journal of Computer Vision 121(2) (2017), 183–208.

43.

Shahdoosti

H.R.

and Hazavei

S.M.

, Combined ripplet and total variation image denoising methods using twin support vector machines[J], Multimedia Tools and Applications 77(6) (2018), 7013–7031.

44.

Mayilvaganan

and Kalpanadevi

, The Survey –Predicting an Education Performance of Students Based on Data Mining Techniques, International Journal of Innovative Research in Computer and Communication Engineering (IJIRCCE) 5(6), 2017.

45.

Wang

, Peng

, Zhao

, et al., Hyperspectral image restoration via total variation regularized low-rank tensor decomposition[J], IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 11(4) (2018), 1227–1243.

46.

Wan

, Gu

, Qian

, et al., Total Variation Regularization Term-Based Low-Rank and Sparse Matrix Representation Model for Infrared Moving Target Tracking[J], Remote Sensing 10(4) (2018), 510.

47.

Gilboa

and Osher

, Nonlocal linear image regularization and supervised segmentation[J], Multiscale Modeling & Simulation 6(2) (2007), 595–630.

48.

Khalid

, Sethu

P.S.

and Sethunadh

, Optimised surfacelet transform based approach for video denoising[C], Inventive Computation Technologies (ICICT), International Conference on, IEEE 3 (2016), 1v5.

49.

Khalid

, Sethu

P.S.

and Sethunadh

, Video denoising using surfacelet transform[C], Recent Trends in Electronics, Information & Communication Technology (RTEICT), IEEE International Conference on IEEE (2016), 1502–1505.

50.

Candes

E.J.

, Wakin

M.B.

and Boyd

S.P.

, Enhancing sparsity by reweighted l1 minimization, Journal of Fourier Analysis and Applications 14(5-6) (2008), 877–905.

51.

Zhang

, Liu

, Xiong

, Ma

and Zhao

, Improved total variation based image compressive sensing recovery by nonlocal regularization, in Proc IEEE Int Symp Circuits Syst (2013), pp. 2836–2839.

52.

Zhang

, Zhao

, Xiong

, Ma

and Gao

, Image compressive sensing recovery via collaborative sparsity, IEEE J Emerg Sel Topics Circuits Syst 2(3) (2012), 380–391.

53.

Mun

and Fowler

J.E.

, Block compressed sensing of images using directional transforms, in Proc 16th IEEE Int Conf Image Process (ICIP) (2009), pp. 3021–3024.

54.

Rajasekar

and Karthikeyan

Dr. , Accessing Secured Public Data Storage Using Secure Index in Cloud Computing Environment, International Journal of Innovative Research in Computer and Communication Engineering (IJIRCCE) 5(7), 2017.

55.

Chen

, Tramel

E.W.

and Fowler

J.E.

, Compressed-sensing recovery of images and video using multihypothesis predictions, in Proc 45th Asilomar Conf Signals Syst Comput (2011), pp. 1193–1198.

56.

Dong

, Shi

, Li

, Ma

and Huang

, Compressive sensing via nonlocal low-rank regularization, IEEE Trans on Image Processing (TIP) 23(8) (2014), 3618–3612.

57.

, Win

, Jing

and Zhang

, An efficient augmented Lagrangian method with applications to total variation minimization, Comput Optim Appl 56(3) (2013), 507–530.

58.

Dong

, Shi

, Ma

and Li

, Image Restoration via Simultaneous Sparse Coding: Where Structured Sparsity Meets Gaussian Scale Mixture,pp, International Journal of Computer Vision (IJCV) 114(2) (2015), 217–232.

59.

, Chen

, Lin

, et al., Nonlocal Joint Regularizations Framework with Application to Image Denoising [J], Circuits Systems and Signal Processing 35 (2016), 2932–3942.

60.

Zuo

, Zhang

, Song

and Zhang

, Texture Enhanced Image Denoising via Gradient Histogram Preservation.