Combination of anisotropic diffusion and non-subsampled shearlet transform for image denoising

Abstract

In this work, a non-subsampled shearlet transform (NSST) based anisotropic diffusion method is proposed. In the proposed method, the NSST transform is firstly applied to the noisy image to provide several scale and directional components. Then, the NSST coefficients are classified into textured regions and noise-related ones by using Sparse Un-mixing by the variable Splitting and Augmented Lagrangian (SUnSAL) classifier. Subsequently, an energy function is formed by the noise-related coefficients to be minimized by diffusion equations. The noisy free image is approximated from the denoised coefficients obtained by the anisotropic diffusion method and textured coefficients which are remained unchanged. Visual and quantitative assessments demonstrate that the proposed method outperforms the state-of-art denoising methods in terms of noise removal and detail preservation.

Keywords

Anisotropic diffusion image denoising edge recovering non-subsampled shearlet transform classification

1 Introduction

A digital image is inevitably corrupted by noise in the process of transmission, coding and receiving. Image denoising methods aim at recovering the unknown initial image from a contaminated observation or noisy measurement. Up to now, numerous image denoising methods including, Gaussian smoothing filtering [1], anisotropic diffusion [2], bilateral filtering [3], total variation minimization [4] and multiresolution wavelet/contourlet/shearlet thresholding method [5 –7] have been explored for recovering the “true” image and removing the noise. These methods adopt different strategies to use the local information of the noisy image for removing the noise.

Nowadays, multiresolution approaches lead to convenient tools in signal processing applications such as, data compression and signal denoising. To find an efficient representation of objects with discontinuities along curves, several researchers have introduced different multiresolution approaches namely, the wavelet transform [8], contourlet transform [9] and shearlet transform [10]. But these transforms are not shift-invariant, which can be a problem in signal analysis, pattern recognition or, as in our case, image denoising [11]. The main drawback of these transforms is that additional edges or structures called artifacts are produced by these methods. To tackle this problem, the shift-invariant version of these transforms i.e., non-subsampled wavelet transform, non-subsampled contourlet transform (NSCT) and non-subsampled shearlet transform (NSST) have been proposed in Refs. [11 –13].

The wavelet method is suitable for one dimensional objects, because this transform does not provide directional components. Among the multiresolution approaches providing directional components, NSST has several advantages. It has more sparse coefficients and a lower computational complexity compared with NSCT [14]. There is no limitation on the number of directions in the NSST to which images are decomposed. In addition, unlike the NSCT, NSST guarantees sparse approximations for piecewise smooth images [15]. Therefore, the NSST has been used to decompose the noisy image in this paper.

Since white noise is distributed evenly over all the wavelet/contourlet/shearlet coefficients, removing the coefficients with small magnitude decreases most of the noise energy while preserving most of the image energy. So, it is understandable why as simple an operation as thresholding in the wavelet/contourlet/shearlet domain can reduce the effect of the noise while preserving image information [16].

Because of its simplicity, the thresholding method is widely used for image denoising in the multiresolution based methods. But the thresholding method presents oscillations in the vicinity of the image’s discontinuity which is similar to the Gibbs phenomena exhibited by the Fourier thresholding [17]. Because of this similarity, they are called pseudo-Gibbs phenomena. Although using the soft thresholding method instead of hard one enables us to partially reduce the pseudo-Gibbs phenomena, this method lowers all of the coefficients and consequently the local average is not preserved leading edges to be eroded [17].

Anisotropic diffusion method follows a completely different approach. To keep edges and textures, the anisotropic diffusion method can be performed using I_t = ∇ . (c (x, y, t) ∇ I), where c (x, y, t) = g (|| ∇ I (x, y, t) ||), in which g is a monotonically decreasing function. As a consequence, for the textured regions with high gradient pixels, c (x, y, t) is small and therefore gets less diffused. On the contrary, for smooth regions with low gradient pixels, c (x, y, t) has a higher value, and these pixels get blurred with neighboring pixels [18]. So, the anisotropic diffusion method smoothes out noise while preserving the informative features.

In some recent papers, several researchers use learning methods e.g., support vector machines (SVMs) to classify the coefficients provided by the multiresolution transforms, into two classes, i.e. clean and noisy [6 , 19]. After determining the noisy coefficients, they are denoised by the soft thresholding approach. Then, the inverse transform is applied to the clean coefficients and the denoised version of the noisy coefficients to get the noise free image. The main advantage of the SVM classifier over other conventional classifiers is that it is very successful when there are a limited number of training samples. To cope with the challenging task of classification in the small sample size situation, several classifiers have been proposed to improve the SVM based classifiers. The SUnSAL classifier is one of the state of the art classifiers aiming at representing most new observations with linear combinations of atoms (training samples), which are selected from an over-complete training dictionary [20]. It has been experimentally shown that, the SUnSAL classifier is more successful than the SVMs based classifiers even in the small sample size satiation. Therefore, we have made use of this classifier in our methodology.

The denoising algorithm presented in this paper combines the NSST and the anisotropic diffusion method by a new strategy. Although, a hybrid method combining the shearlet and anisotropic diffusion, has been proposed for denoising in Ref. [21], we improve it in this paper by classifying the NSST coefficients and forming a new cost function to be minimized. The proposed method takes advantage of multiresolution approach while the artifacts have been reduced because of the anisotropic diffusion used in our method.

This article is organized into five sections.Section 2 reviews the NSST, the SUnSAL classifier and the anisotropic diffusion method. In Section 3, the proposed method is presented. The experimental results and discussions are described in Section 4. Finally, Section 5 is devoted to the conclusions.

2 A review of methods

This section reviews the NSST, the SUnSAL classifier and the anisotropic diffusion method as presented in Refs. [13 , 22].

2.1 NSST

The NSST is a shift-invariant, multi-scale and multi-directional expansion demonstrating low computational complexity compared with NSCT. Moreover, unlike other multiresolution transforms, NSST transform can be decomposed into any number of directions. Therefore, using NSST in the image denoising application provides more sparse coefficients by which the denoising algorithm has a higher performance. Before describing NSST, the discrete shearlet transform should be explained. So, in this subsection the discrete shearlet transform and NSST are reviewed.

A) Discrete shearlet transform

When dimension n is 2, the affine systems with composite dilations are defined as following:

$\begin{matrix} A_{AS} (Ψ) = {Ψ_{j, l, k} (x) \\ = | \det A |^{j / 2} Ψ (S^{l} A^{j} x - k) : j, l \in Z, k \in Z^{2}} \end{matrix}$ (1) in which Ψ ∈ L² (R²), S and A are 2×2 invertible matrices where, |detS|=1. For the appropriate choices of ψ, A_AS (Ψ) forms a tight frame for L² (R²) and the elements of this system are called composite wavelet, i.e.: $\sum_{j, l, k} | < f, Ψ_{j, l, k} > |^{2} = | | f | |^{2}$ (2) where f is a function on which the shearlet is implemented. If A is the parabolic scaling matrix and S is the shearing matrix, for each a > 0 (a is 2^-j) and s ∈ R (s = - l) we have: $A = (\begin{matrix} a & 0 \\ 0 & \sqrt{a} \end{matrix}), S = (\begin{matrix} 1 & s \\ 0 & 1 \end{matrix})$ (3)

The first matrix, i.e. A is an anisotropic dilation matrix which controls the ‘scale’ of the shearlets and guarantees that the frequency support of the shearlets becomes increasingly elongated at finer scales, while the second matrix, i.e. S is a shearing matrix controlling the orientation of the shearlets. It is usually assumed a = 4 and s = 1, such that $A = A_{0} = (\begin{matrix} 4 & 0 \\ 0 & 2 \end{matrix})$ and $S = S_{0} = (\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix})$ .

For $\forall ξ = (ξ_{1}, ξ_{2}) \in {\hat{R}}^{2}, ξ_{1} \neq 0,$ assume ${\hat{Ψ}}^{0} (ξ)$ is: ${\hat{Ψ}}^{0} (ξ) = {\hat{Ψ}}^{0} (ξ_{1}, ξ_{2}) = {\hat{Ψ}}_{1} (ξ_{1}) {\hat{Ψ}}_{2} (\frac{ξ_{2}}{ξ_{1}})$ (4) in which ${\hat{Ψ}}_{1} \in C^{\infty} (R)$ is a wavelet and $supp {\hat{Ψ}}_{1} \subset [- 1 / 2, - 1 / 16] \cup [1 / 16, 1 / 2]$ . In addition, ${\hat{Ψ}}_{2} \in C^{\infty} (R)$ , and $supp {\hat{Ψ}}_{2} \subset [- 1, 1]$ . So, one may conclude ${\hat{Ψ}}^{0} \in C^{\infty} (R)$ and $supp {\hat{Ψ}}^{0} \subset [- 1 / 2, 1 / 2]^{2}$ .

For any (ξ₁, ξ₂) ∈ C₀, where $C_{0} = {(ξ_{1}, ξ_{2}) \in {\hat{R}}^{2} : | ξ_{1} | \geq 1 / 8, | ξ_{2} | \leq 1}$ called horizontal cone, the functions ${{\hat{ψ}}^{0} ({ξ A}_{0}^{- j} S_{0}^{- l})}$ form a tiling of C₀ (see Fig. 1). This property implies that the collection $\begin{matrix} {Ψ_{j, l, k}^{0} {(x) = 2}^{\frac{3 j}{2}} Ψ^{0} (S_{0}^{l} A_{0}^{j} x - k) : \\ j \geq 0, - 2^{j} \leq l \leq - 2^{j} - 1, k \in Z^{2}} \end{matrix}$ (5) forms a Parseval frame for $L^{2} {(C}_{0})^{\lor} = {f \in L^{2} {(R}^{2}) : supp \hat{f} \subset C_{0}}$ . Analogously to C₀, a Parseval frame for L² (C₁) ^∨ can be constructed, where C₁ is the vertical cone: $C_{1} = {(ξ_{1}, ξ_{2}) \in {\hat{R}}^{2} : | ξ_{2} | \geq 1 / 8, | \frac{ξ_{1}}{ξ_{2}} | \leq 1}$ (6)

Suppose $A_{1} = (\begin{matrix} 2 & 0 \\ 0 & 4 \end{matrix})$ , $S_{1} = (\begin{matrix} 1 & 0 \\ 1 & 1 \end{matrix})$ and ${\hat{Ψ}}^{1} (ξ) = {\hat{Ψ}}^{1} (ξ_{1}, ξ_{2}) = {\hat{Ψ}}_{1} (ξ_{1}) {\hat{Ψ}}_{2} (\frac{ξ_{1}}{ξ_{2}})$ . Then, the Parseval frame for L² (C₁) ^∨ is as follows, $\begin{matrix} {Ψ_{j, l, k}^{1} {(x) = 2}^{\frac{3 j}{2}} Ψ^{1} (S_{1}^{l} A_{1}^{j} x - k) : \\ j \geq 0, - 2^{j} \leq l \leq - 2^{j} - 1, k \in Z^{2}} \end{matrix}$ (7)

B) Non-subsampled shearlet transform (NSST)

The shift invariant version of the shearlet transform was proposed by Easley in 2008 [13]. The NSST eliminates the down-sampler and up-sampler operators in the structure of the shearlet transform. This transform uses the non-subsamples Laplacian pyramid (NSLP) with several combinations of shearing filters. The process by which the scales are produced is as follows: ${Scale}_{j + 1} = h_{j}^{s} I = (h_{j}^{l} \prod_{k = 1}^{j - 1} h_{k}^{h}) I$ (8) where I is a function (I is an image in the image denoising application), Scale_j+1 is the subband at scale j + 1, and $h_{k}^{l}$ and $h_{j}^{h}$ denote low-pass and high-pass filters of NSLP corresponding to the scales jand k.

Given an N × N image $I_{a}^{0}$ and the number of direction D, the steps of NSST in the scale j are described as follows:

Use the low-pass and high-pass filters of NSLP to decompose $I_{a}^{j - 1}$ into a low-pass image $I_{a}^{j}$ and a high-pass image $I_{d}^{j}$ . Note that removing the downsampler results in $I_{a}^{j}$ and $I_{d}^{j}$ with size N × N.

Perform FFT to $I_{d}^{j}$ to obtain ${\hat{I}}_{d}^{j}$ in pseudo polar grid, then compute ${PI}_{d}^{j}$ .

To get ${{\hat{I}}_{d, k}^{j}}_{k = 1}^{D}$ , apply a band-pass filtering to ${PI}_{d}^{j}$ .

Performing inverse FFT to ${{\hat{I}}_{d, k}^{j}}_{k = 1}^{D}$ , gives NSST coefficients ${I_{d, k}^{j}}_{k = 1}^{D}$ in pseudo polar grid.

Figure 2 shows the results of applying the two level NSST on the Zoneplate image, yielding a low-pass subband (Fig. 2(b)) and several highpass subbands. In this figure, the numbers of shearing directions are chosen to be 4 and 8 from coarser to finer scale.

2.2 SUnSAL classifier

The sparse representation is based on the idea that all the test samples can be well-approximated by a (sparse) linear combination of atoms selected from an overcomplete training dictionary [20]. Given a training dictionary, denoted as A = {x₁, …, x_n} ∈ R^n×m, in which x_i ∈ R^1×m is a training sample, n is the number of training samples and m denotes the dimension of each training sample (i.e., the number of features), we want to approximate each new observation by this dictionary. Let the dictionary contain a total of c different classes. So, the dictionary can be organized as A = [A₁, …, A_c], where A_k = {x_{k
₁}, …, x_{k
_{n
_k}}} and A_k includes the training samples of class k in its column, n_k is the number of training samples in class k, and $\sum_{k = 1}^{c} n_{k} = n$ . In supervised classification of data, a conventional scenario is that a limited set of labeled training samples for each class exists, and we use this information to train a classifier, which is then used to classify unlabeled samples. Given a new observation (an unlabeled sample) denoted as x_i, sparse representation aims at representing it by linear combination of the training samples (atoms) in the dictionary A as follows:

$\begin{matrix} x_{i} \approx x_{1} a_{1} + x_{2} a_{2} + \dots + x_{n} a_{n} \\ = [x_{1} x_{2} \dots x_{n}] [a_{1} a_{2} \dots a_{n}]^{T} \\ = [A_{1}, \dots, A_{c}] [a_{1}^{T}, \dots, a_{c}^{T}]^{T} = A a + ɛ \end{matrix}$ (9) in which $a = [a_{1}^{T}, \dots, a_{c}^{T}]^{T} \in R^{n \times 1}$ is a sparse vector (i.e., most of its elements are zero), a _i ∈ R^n_i×1 is the vector of regression coefficients associated with class i, and ɛ denotes the error of representation. It is supposed that x_i can be well approximated by A_i a _i, such that a _j = 0, for j ≠ i. In summary (1) For a given unlabeled sample x_i, a should satisfy x_i = A a + ɛ. (2) Vector a should be as sparse as possible. So, solving the following optimization problem gives the sparse vector a for which we are looking: $\hat{a} = \arg min | | a {| |}_{0} subject to x_{i} = E e$ (10) where || a ||₀ denotes the l₀-norm in which the nonzero elements are counted. Owing to the fact that the noise exists, Equation (10) is replaced by: $\hat{a} = \arg \min | | a {| |}_{0} subject to | | x_{i} - A a {| |}_{2} \leq δ$ (11) where δ presents an error tolerance. Equation (11) is very difficult to be solved due to the presence of l₀-norm in this equation. To cope with this problem, recently proposed algorithms such as basis pursuit (BP) and orthogonal matching pursuit (OMP) are usually used. To solve the problem, BP substitutes the l₀-norm with the l₁-norm, while OMP solves Equation (11) by a greedy strategy in the several iterations [23, 24]. Substituting l₀-norm with the l₁-norm gives: $\hat{a} = \arg \min | | x_{i} - A a {| |}_{2}^{2} + τ_{1} | | a | |_{1}$ (12)

Note that Equation (12) is a convex problem. Therefore, it can be solved using linear programming solvers. The Lagrange multiplier denoted by τ₁ creates a tradeoff between the remained error and sparse solution. Analogously to the results of Ref. [20], we have experimentally observed that the nonnegativity constraint (NC) i.e., a ≥ 0 improves the PSNR of the approximated image. Due to the fact that there is a positive dependency/correlation among the adjacent pixels of the non-subsampled transforms, representation of training samples should yield positive values for coefficients a . Therefore, the optimal sparse vector can be obtained as follows: $\hat{a} = \arg \min | | x_{i} - A a {| |}_{2}^{2} + τ_{1} | | a | |_{1}, a > 0$ (13) which is an l₂ - l₁ optimization problem. The sparse unmixing by variable splitting and augmented Lagragian (SUnSAL) method is an appropriate tool for solving Equation (13). As it is expected, due to the noise involved in the problem, the solution of Equation (13) is not ideal. Therefore, to assign the class label of the sample x_i, the following residual-based criterion is used: ${cl \hat{a} ss (x}_{i}) = \arg \min_{j \in {1, \dots, c}} | | x_{i} - A_{j} a_{j} {| |}_{2}$ (14)

2.3 Anisotropic diffusion

The idea of image denoising by the nonlinear diffusion method can be explained as follows [22]. Let I be the noisy image. This image is known to be the sum of the noise free image and some Gaussian noise n (x, y): $I (x, y) = f (x, y) + n (x, y), (x, y) \in Ω .$ (15) where, Ω ⊂ R² denotes a rectangle. The diffusion process is as follows: $\frac{\partial I}{\partial t} = \nabla \cdot (g | \nabla I | \nabla I)$ (16) with the initial condition: $I (x, y, 0) = I (x, y)$ (17) in which the time t is the scale parameter for filtering. The diffusion function g (x) = const corresponds to linear diffusion leading to a strong smoothing of I with increasing t. Typically, g (x) is selected to be a non-negative decreasing function with $lim_{x \to \infty} g (x) = 0$ . The diffusivity g (x) controls the smoothing process by admitting strong diffusion if the gradient ∇I is small which is likely to correspond to the location of noise, and by slowing down (or even stop) the smoothing for large gradients which are likely to correspond to the location of an edge.

The most well-known diffusivity, is Perona–Malik diffusivity g (x) =1/(1 + x²/γ²) which is used in the proposed method.

3 Proposed method

The goal of this section is to propose a model combining NSST and nonlinear diffusion for discontinuity-preserving denoising. Let I (x, y) be an observed discrete noisy image on the square Ω ⊂ R² known to be the sum of the original image and noise. We first apply NSST to this noisy image. In the conventional hybrid methods, such as curvelet–diffusion [22] and shearlet-diffusion [21], the coefficients of detail subbands of the transform are divided into two sets by a threshold value. In these methods, coefficients with values more than the threshold are preserved while the coefficients with values lower than the threshold are denoised by the diffusion filtering. This idea preserves the coefficients with higher values which are likely to correspond to the location of textured regions. The reason is that these coefficients contain informative features of the image and denoising these coefficients can erode the edges of the image. However, determining the coefficients belonging to the textured regions cannot be done with only a simple thresholding operation. So, to improve the quality of the denoised image, we adopt a different strategy in which the coefficients belonging to the textured regions and noise-related ones, are determined by the SUnSAL classifier. Then, the coefficients of the second class, i.e., the noise-related class, are denoised by the diffusion filtering.

In the NSST, pixels of the original image spatially relate to the coefficients in the transform domain. The spatial relationships in the NSST domain represent the edges and textures of the image which should be preserved in the denoising process. A natural image is formed by regular features and performing amultiresolution transform such as NSST on the image results in sparse, spatially joined coefficients, representing edges and textures of the original image. The spatial dependency/contiguity of the multiresolution transform coefficients is known as spatial regularity [7]. Some researchers have used the similar concept (i.e., signal regularity) instead of spatial regularity for determining the important coefficients in their proposed method [19]. Although the spatial dependency of the features is represented by statistical models such as Markov random sequences in some methods [25], these forms of representation have two defects. First, the computational complexity of these methods is high. Second, their models cannot describe the geometry of the features truly. To remove these defects, the spatial dependency of the features is described by the connectivity of NSST coefficients in our strategy. Due to the fact that the NSST realizes that spatial dependency between its coefficients, the subbands of the NSST are generally formed by the connected coefficients. This concept is utilized in our methodology to determine the textured regions and noise-related ones. After determining the noise-related coefficients, they are denoised by the diffusion filtering in our method to prevent the artifacts usually produced by the soft thresholding method. The steps of the proposed hybrid method are asfollowing:

(1) Perform the NSST on the noisy image I to obtain L levels and O orientations. One approximate NSST subband denoted by W⁰ and several detail subbands denoted by W_i,j (i = 1, 2, …, L, j = 1, 2, …, O) are provided here.

(2) Select the training samples to be used as atoms, and find the feature vector of these samples for detail subbands of the NSST i.e., W_i,j. For this purpose, firstly obtain the preliminary binary map B_i,j [x, y]. Then compute the support value S_i,j [x, y] from the binary map, in which x and y determine the position of coefficient. The binary map B_i,j [x, y] is obtained by the following formula: $B_{i, j} [x, y] = {\begin{matrix} 1, when | W_{i, j} [x, y] | > T \\ 0, else \end{matrix}$ (18) where threshold T determines the coefficients with high values which are likely to correspond to the location of an edge (see Fig. 3(b)). For each coefficient, sum the values of the neighboring coefficients of the binary map to obtain the support value S_i,j [x, y]. In fact, the support value determines the number of neighboring coefficients with values larger than the threshold T. So, one expects that coefficients with high support values exhibit features of the image and coefficients with low support values represent noise. This step has been shown in Fig. 3(c). Subsequently, k coefficients with the max support value are selected as training samples for textured coefficients, and k coefficients with support value zero are randomly selected as training samples for noise-relatedcoefficients. After selecting the training samples, the class number is assigned to them which is 1 for textured coefficients and 0 for noise-related coefficients.

(3) Use the SUnSAL classifier to classify all the coefficients of detail subbands W_i,j (i = 1, 2, …, L, j = 1, 2, … , O), where the training samples are used as atoms.

(4) Now, the following idea is used. We want to keep the low frequency coefficients and the important NSST coefficients (labeled as textured regions) almost untouched, but we want to denoise the NSST coefficients which have been labeled as noisy coefficients, by the diffusion procedure in such a way, that the image is smoothed. Therefore, we will not apply the nonlinear diffusion process directly to I, but only to the difference image I_d = I - I₀, where I₀ is obtained by applying the inverse NSST to the set of coefficients formed by low frequency coefficients and textured coefficients. After some iteration steps of the diffusion scheme, the smoothed difference image ${\hat{I}}_{d}$ is added to I₀ in order to obtain the final result. The proposed model can be formulated as follows: $\frac{\partial I_{d}}{\partial t} = \nabla \cdot (g | \nabla I_{d} | \nabla I_{d})$ (19) with I_d (x, y, 0) = I - I₀ denoting the difference image.

4 Experimental results and discussions

We have conducted some experiments to assess the performance of the proposed method quantitatively and visually. In the quantitative comparison, noisy images were simulated by adding white Gaussian noise with various standard deviations (20, 30, 40, 50, 60, 70) to some standard test images namely, Lena, Boat, Barbara, Cameraman, and a set of ordinary images from the Berkeley dataset including Castle, Man and Children (see Fig. 4). The noise was then removed from these images using several algorithms and the PSNR results were calculated. Four-level nonsubsampled shearlet transform (NSST) with 2, 2, 3, 3 directions, respectively, was used to decompose the noisy images.

In the first step, we should determine the free parameters used in our algorithm. Our experiments showed that the optimal regularization parameter is τ₁ = 3 ×10^-5. In addition, the number of training samples k (see step 2 of the proposed algorithm) is set to be 24, because the improvement is negligible for the larger number of k. In the next step, we should determine the threshold value T, to be used in the proposed method. For this purpose, the idea called oracle [19] is usually used in the image processing field, in which several images are corrupted by noise and subsequently, the noise is removed by a series of T and the optimal threshold T maximizing the PSNR value is obtained to be used for noise reduction in general. Here, five well-known training images, namely, Stream and bridge, Fruits, Goldhill, Couple and Baboon have been used. The obtained optimal threshold T is 11. In addition the time step used to iteratively solve Equation (19), was set to be 0.05.

In the implementations, we have used the publicly codes available for BLS-GSM [26], ProbShrink [27] and NL-means [28], but quantitative results of ShearletSVM have been reported from Ref. [7] just in the table.

Visual inspections are necessary to evaluate the performance of image denoising algorithms, because artificial edges cannot be well-quantified by the PSNR value. The results of the denoising algorithms for the Lena data have been represented in Fig. 5. Zoomed areas of Fig. 5 have also been shown in Fig. 6. In addition, the visual performances of the denoised images for the Cameraman data have been represented in Fig. 7. The standard deviation of noise is 20 in this experiment. From these figures (specifically Fig. 6(c)), one can easily observe that, in the BLS-GSM method, the details and edges have been eroded, but the noise has been reduced satisfactorily. As it is obvious, the result of the ProbShrink method has preserved the edges better than that of the BLS-GSM method, but the noise is still remained in the result of this method. In addition, a considerable amount of artifacts has been produced by this method (see Figs. 5, 6, 7(d)). As can be seen from these figures, the result of NL-means method seems pretty noisy, such that the level of the noise in this method is higher than that in the ProbShrink method. Moreover, the edges and textures have been blurred in this method. One can observe that the proposed denoising method preserves the edges and textures better than other methods, while removing the noise almost completely (see Figs. 5, 6, 7(f)).

Finally, the performance of these methods was assessed quantitatively. The values of PSNR index for the proposed method and other mentioned methods have been reported in Table 1. This table demonstrates that the proposed method outperforms other powerful denoising methods quantitatively.

5 Conclusion

A new image denoising method was proposed combining NSST and anisotropic diffusion. In order to preserve features of the image, after applying the NSST to the noisy image, the obtained coefficients are classified into two classes i.e., edge-related coefficients and noise-related ones by using the SUnSAL classifier, and only the coefficients of the second class are used in the denoising procedure. To prevent artifacts, we denoised the noise-related coefficients by the anisotropic diffusion method. Some standard grayscale images and several ordinary images selected from the Berkeley dataset were denoised by the proposed method and five powerful algorithms. The comparison of the denoising results obtained with the proposed method, and with the several powerful methods, demonstrates the efficiency of our new image denoising approach. The visual quality of denoised images obtained by the proposed method is moreover characterized by fewer artifacts than the other methods.

References

Lindenbaum

, Fischer

and Bruckstein

A.M.

, On Gabor’s contribution to image enhancement, Pattern Recognit27 (1994), 1–8.

Perona

and Malik

, Scale-space and edge detection using anisotropic diffusion, IEEE Trans Pattern Anal Mach Intell12 (1990), 629–639.

Zhang

and Gunturk

B.K.

, Multiresolution bilateral filtering for image denoising, IEEE Transactions on Image Processing17 (2008), 2324–2333.

Oliveira

J.P.

, Bioucas-Dias

J.M.

and Figueiredo

M.A.T.

, Adaptive total variation image deblurring: A majorization–minimization approach, Signal Process89 (2009), 1683–1693.

Donoho

D.L.

, De-noising by soft-thresholding, IEEE Trans Inform Theory41 (1995), 613–627.

Shahdoosti

H.R.

and Khayat

, Image denoising using sparse representation classification and non-subsampled shearlet transform, Signal, Image and Video Processing, DOI: 10.1007/s11760-016-0862-0

Yang

H.Y.

, Wang

X.Y.

, Niu

P.P.

and Liu

Y.C.

, Image denoising using nonsubsampled shearlet transform and twin support vector machines, Neural Networks57 (2014), 152–165.

Mallat

S.G.

, A theory for multiresolution signal decomposition: The wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence11 (1989), 674–693.

M.N.

and Vetterli

, The contourlet transform: An efficient directional multiresolution image representation, IEEE Transactions on Image Processing14 (2005), 2091–2106.

10.

Gomathi

and Kumar

, Inpainting for satellite imagery using thin plate spline radial basis function neural networks in shearlet domain, Journal of Intelligent & Fuzzy Systems27 (2014), 2391–2398.

11.

Nunez

, Otazu

, Fors

, Prades

, Pala

and Arbiol

, Multiresolution-based image fusion with additive wavelet decomposition, IEEE Transactions on Geoscience and Remote Sensing37 (1999), 1204–1211.

12.

Da Cunha

A.L.

, Zhou

and Do

M.N.

, The nonsubsampled contourlet transform: Theory, design, and applications, IEEE Transactions on Image Processing15 (2006), 3089–3101.

13.

Easley

, Labate

and Lim

W.Q.

, Sparse directional image representations using the discrete shearlet transform, Applied and Computational Harmonic Analysis25 (2008), 25–46.

14.

Guorong

, Luping

and Dongzhu

, Multi-focus image fusion based on non-subsampled shearlet transform, IET Image Processing7 (2013), 633–639.

15.

Lim

W.Q.

, Nonseparable shearlet transform, IEEE transactions on image processing22 (2013), 2056–2065.

16.

Wang

X.Y.

, Yang

H.Y.

and Fu

Z.K.

, A New Wavelet-based image denoising using undecimated discrete wavelet transform and least squares support vector machine, Expert Systems with Applications37 (2010), 7040–7049.

17.

Durand

and Froment

, Reconstruction of wavelet coefficients using total variation minimization, SIAM J Sci Comput24 (2003), 1754–1767.

18.

Jianwei

and Plonka

, Combined curvelet shrinkage and nonlinear anisotropic diffusion, IEEE Transactions on Image Processing16 (2007), 2198–2206.

19.

Balster

E.J.

, Zheng

Y.F.

and Ewing

R.L.

, Feature-based wavelet shrinkage algorithm for image denoising, IEEE Transactions on Image Processing14 (2005), 2024–2039.

20.

Song

, Li

, Dalla Mura

, Peijun

, Plaza

, Bioucas-Dias

J.M.

, Benediktsson

J.A.

and Chanussot

, Remotely sensed image classification using sparse representations of morphological attribute profiles, IEEE Transactions on Geoscience and Remote Sensing52 (2014), 5122–5136.

21.

Easley

, Labate

R.D.

and Colonna

, Shearlet-based total variation diffusion for denoising, IEEE Transactions on Image Processing18 (2009), 260–268.

22.

and Plonka

, Combined curvelet shrinkage and nonlinear anisotropic diffusion, IEEE Transactions on Image Processing16 (2007), 2198–2206.

23.

Chen

, Donoho

and Saunders

, Atomic decomposition by basis pursuit, SIAM J Sci Comput43 (2001), 129–159.

24.

Tropp

and Gilbert

, Signal recovery from random measurements via orthogonal matching pursuit, IEEE Trans Inf Theory53 (2007), 4655–4666.

25.

Qiu

and Mukherjee

P.S.

, Edge structure preserving image denoising, Signal Processing90 (2010), 2851–2862.

26.

Portilla

, Strela

, Wainwright

M.J.

and Simoncelli

E.P.

, Image denoising using scale mixtures of gaussians in the wavelet domain, IEEE Trans Image Process12 (2003), 1338–1351.

27.

Pižurica

and Philips

, Estimating the probability of the presence of a signal of interest in multiresolution single- and multiband image denoising, IEEE Trans Image Process15 (2006), 645–665.

28.

Buades

, Coll

and Morel

J.M.

, A non-local algorithm for image denoising, In Computer Vision and Pattern Recognition, CVPR 2005, IEEE Computer Society Conference on2 (2005), 60–65.