Multimed Tools Appl DOI 10.1007/s11042-015-3154-8

An adaptive edge-preserving image denoising technique using patch-based weighted-SVD filtering in wavelet domain Paras Jain 1 & Vipin Tyagi 1

Received: 21 May 2015 / Revised: 7 December 2015 / Accepted: 9 December 2015 # Springer Science+Business Media New York 2015

Abstract Image denoising has always been one of the standard problems in image processing and computer vision. It is always recommendable for a denoising method to preserve important image features, such as edges, corners, etc., during its execution. Image denoising methods based on wavelet transforms have been shown their excellence in providing an efficient edge-preserving image denoising, because they provide a suitable basis for separating noisy signal from the image signal. This paper presents a novel edge-preserving image denoising technique based on wavelet transforms. The wavelet domain representation of the noisy image is obtained through its multi-level decomposition into wavelet coefficients by applying a discrete wavelet transform. A patch-based weighted-SVD filtering technique is used to effectively reduce noise while preserving important features of the original image. Experimental results, compared to other approaches, demonstrate that the proposed method achieves very impressive gain in denoising performance. Keywords Wavelet transform . Singular value decomposition (SVD) . Edge detection . Adaptive filtering

1 Introduction In modern times, images and videos have become an integral part of our life. Applications have been extended from general documentation of events and visual communication to the more serious surveillance and medical fields. This has increased the demand for accurate and visually impressive images. However, digital images may often get corrupted by noise during the process of acquisition and transmission. This

* Vipin Tyagi [email protected] 1

Jaypee University of Engineering and Technology, Raghogarh, Guna, MP 473226, India

Multimed Tools Appl

form of corruption may degrade the visual quality of the image. The major challenge for a filtering technique is to suppress noise without affecting the important image structures and details. Several approaches [17, 20, 39] for noise reduction have been proposed in last few decades, many of them are based on linear spatial domain filters. Linear spatial filters (for example, Gaussian filters [39]) are simple and easily implementable and usually smooth the data to reduce noise effects; however, they can also result in blurring of important image structures such as edges [17]. Furthermore, most linear filtering techniques proposed so far need some prior knowledge about the noise and the image characteristics. Mostly, such details are not available and may be hard to estimate from the input noisy data. To overcome these problems of linear filters, filtering based on non-linear edgepreserving methods have become the main stream in the field and research is still going on in this direction [21]. These non-linear edge-preserving filtering methods can preserve the important image features while suppressing the undesirable noise. Some of these techniques are based on partial differential equations (PDE’s) and variational models. For example, non-linear/anisotropic diffusion (AD) [31] was designed to overcome blurring issues of Gaussian filter [39] by smoothing the image only in the direction orthogonal to the gradient. The regularization methods based on the total variation (TV) [35, 36] were designed to smooth the homogenous regions of the image but not its edges. In [41], Tomasi and Manduchi proposed a simple, non-iterative and local smoothing filter known as the bilateral filter which smooths images by means of a non-linear combination of nearby image values. The non-local means (NLM) [3] is the first filter that makes use of the self-similarity in the whole image. NLM obtains a denoised patch by weighted averaging all other patches in the same image. It is an extension of the bilateral filter [41] in the sense of replacing the Euclidean distance between two pixels with the weighted Euclidean distance between two patches. While the denoising approaches mentioned so far come under the category of spatial domain denoising methods, a vast section of image denoising literature is devoted to transform domain methods. In [38], Shao et al. introduced a new taxonomy of state-ofthe-art image denoising techniques which helps us to choose an appropriate technique based on image representation. They have provided a good comparative study on spatial domain, transform domain, and dictionary learning based denoising approaches. Inspired by non-local filtering in [3], Dabov et al. [9] proposed a transform-based block-matching and 3D (BM3D) filtering method which has great potential. Recently, the BM4D is proposed as an extension of the BM3D algorithm for denoising of volumetric data [27]. Qiu et al. [11] proposed a new edge-preserving smoothing filter, called Linear-Local SURE (LLSURE) filter based on a local linear model and the principle of Stein’s Unbiased Risk Estimate (SURE). Singular value decomposition (SVD) is also used for image noise reduction [4, 19, 24, 30, 46]. Cross validation techniques are used for thresholding parameter selection [28, 45]. Bayesian procedures combine inference from data with prior information to estimate thresholding parameters [7, 43]. Furthermore, a different subgroup of transform domain methods which include the wavelet-based denoising methods, exploits the decomposition of the data into a wavelet basis and modifies the wavelet coefficients to denoise the data [2, 8, 11, 22, 23, 26, 32, 40]. In particular, Jain and Tyagi [23] have proposed very recently a new technique for noise reduction in wavelet domain. They presented a new locally adaptive patch-based (LAPB) thresholding scheme that relies on the aggregation of multiple threshold estimates of a

Multimed Tools Appl

wavelet coefficient and involves estimation of thresholding parameters for a coefficient in a local neighborhood. Wavelet transforms show localization in both time and frequency. The localized nature of the wavelet transforms results in denoising with edge preservation [26, 40]. The basic idea regarding the noise suppression based on the wavelet transform is to compute the multiscale wavelet decomposition of the corrupted image and to modify (threshold) the obtained wavelet coefficients by comparing them with a predetermined threshold according to a shrinkage rule. Reconstruction from these modified coefficients then produces the desired denoised image. Wavelet shrinkage is an often used approach for denoising due to its simplicity. There are several studies on thresholding the wavelet coefficients [5, 6, 11–13, 37]. A great challenge in the wavelet shrinkage process is to determine a sufficient threshold value. The threshold should neither be too small nor be much larger. A smaller threshold may result in insufficient noise suppression, whereas a larger threshold will result in over shrinkage of the signal that may suppress important features of the image. The commonly used threshold estimation criteria are VisuShrink (nonadaptive) [12], SureShrink [13] (adaptive), BayesShrink [6, 7] (adaptive) and Cross Validation [28, 45] (adaptive). The shrinkage (thresholding) rule defines the applicability of a threshold. Basically, each wavelet coefficient is compared against a threshold. If the coefficient is larger than the threshold, the coefficient is preserved for later use and at times for modification; otherwise it is assigned zero. The main objective of the thresholding procedure is to preserve large coefficients which represent important signal features, while small coefficients which mainly represent the noise content can be thresholded without affecting the significant image features. Some well-known shrinkage rules are ‘keep-or-kill’ hard thresholding and ‘shrink-or-kill’ soft thresholding [11], hyperbolic function [43], firm thresholding [15], garrote thresholding [14], SCAD thresholding [1]. This paper presents a new method for noise reduction in wavelet domain. An adaptive patch-based weighted-SVD filtering is used as a thresholding scheme to effectively reduce Gaussian noise while preserving important features of the original image. Experimental results show that the proposed method, when compared to well-known state-of-the-art denoising methods, is more suitable for various classes of images corrupted by Gaussian noise. The remainder of this paper is organized as follows. Section 2 gives a brief review about some SVD-based techniques for noise suppression; Section 3 elaborates the proposed method in detail; Section 4 provides the comparison of the results obtained through the proposed approach with other state-of-the-art denoising methods; Section 5 summarizes the conclusions.

2 Review of SVD-based noise suppression SVD of a matrix is one of the most valuable tools in the theory of inverse problems. It has been successfully applied to many image restoration problems. Common applications include regularization, signal estimation, noise reduction, signal detection, etc. [42]. The SVD-based methods, in order to discriminate against noise, decompose the column space of given noisy matrix into dominant and subdominant subspaces; one of them is attributed to the noise-free signal, hence named the signal subspace and the other is attributed to the noise hence named the noise subspace. Usually, these two subspaces are

Multimed Tools Appl

orthogonal to each other which implies that the signal is independent to noise. For image noise suppression, the basic approach is to minimize the energy in the noise subspace whereas retaining the energy in the signal subspace. In the theory of SVD [16], a realvalued M × N matrix X can be decomposed uniquely as X ¼ U ΣV t ¼

N X

t

αi ! u i! vi

ð1Þ

i¼1 t v i , and where U and V are orthogonal matrices with column vectors ! u i and ! Σ = diag(α1, α2, … …, αN) is a diagonal matrix with diagonal elements αi, i = 1,….., N. The diagonal elements of Σ can be arranged in a nonincreasing order and are called the singular values of X. For convenience, we denote the singular values of matrix X as αi(X). With SVD, one can determine the effective rank of a matrix conveniently. In theory, the rank of X is the number of its nonzero singular values. Specifically, if the rank of X is r, where r ≤ min(M, N), then

α1 ≥ α2 ≥ … ≥αr > αrþ1 ¼ … ¼ αN ¼ 0

ð2Þ

In practice, under an additive noise model, the perturbed (noisy) observation matrix can be obtained as: Y ¼ X þ n;

ð3Þ

where n is a random noise perturbation matrix of full rank. As can be seen in Eq. (2) that last N-r singular values of r-rank matrix X are equal to zero, but when we add it with a noise matrix n of full rank, then noise will expand the whole column space and last N-r singular values of observation matrix Y can be small, but not necessarily zero. Let β1 ≥ β2 ≥ … ≥ βn − 1 ≥ βN be the singular values of the observation matrix Y arranged in nonincreasing order. It can be observed that singular values βr + 1 … βN are significantly smaller than those corresponding to the signal, i.e., β1 … βr. The dominant subspace is due to the signal, hence referred to as signal subspace, while the other is referred to as noise subspace. The effective rank of Y will be r if [24] β r > λ > β rþ1

ð4Þ

where 1 ≤ r ≤ min(M, N) and λ = ∥n∥2 is the 2-norm of n [25]. To determine the effective rank of a perturbed matrix Y, Konstantinides and Yao [42] derived tighter threshold bounds for perturbation due to i.i.d. random noise of zero mean and variance σ2. pffiffiffiffiffiffiffiffiffi pffiffiffi cσ≤λ ≤ M N σ;

ð5Þ

where parameter c is estimated from the statistics of signal and noise. Simulation in [25] demonstrated that the simple threshold λ = 3 σ performs most stably under a variety of noise levels. Let Y ¼ U Y Σ Y V tY ¼

XN

β i¼1 i

t ! u Yi! v Yi

ð6Þ

Multimed Tools Appl

be the SVD of noisy matrix Y = X + n. The estimation Xbof the original image X from the noisy observation Y can be computed as [24]: t Xb ¼ U Y Σ Xb V Y ð7Þ  0  0 0 where ∑Xb ¼ diag α1 ; α2 ; ……; αN is a diagonal matrix having the singular values   0 α′1, α′2, … …, α′N as the diagonal elements and αi Xb is the i-th singular value of the denoised matrix Xb. Each singular value α'i is obtained from the singular value βi of the perturbed matrix Y using the hard thresholding function expressed as [24]:  0 βi ; if β i > λ ð8Þ αi ¼ T h ðβi ; λÞ ¼ 0; if β i ≤ λ With the value of λ chosen as per the condition in (4), the estimation Xb in (7) is the rank-r approximation of the rank-N perturbed matrix Y, (r < N), obtained through the concept of truncated SVD in (8). Therefore, the expression for Xb can be rewritten as: Xb ¼

Xr

β i¼1 i

t ! u Yi! v Yi

ð9Þ

The most impressive fact of Eq. (8) is that it allows us to robustly approximate original signal when observation is corrupted with noise. The important feature of truncated SVD for suppressing the noise in an image is to retain the energy in the signal subspace and discard the energy in the noise subspace. The small singular values which usually are the part of the noise subspace, are truncated to zero to obtain the filtered matrix with reduced rank. The crucial point is to determine the threshold which differentiates the signal subspace from the noise subspace. The traditional hard thresholding operation for SVD-based filtering which is expressed in Eq. (8) exhibits some discontinuities and may be unstable or more sensitive to small changes in the data. As an alternative, Cai et al. [4] suggested to use soft thresholding function to determine the truncated singular values, expressed as: 0

αi ¼ T s ðβ i ; λÞ ¼ sgnðβ i Þmaxðjβ i j−λ; 0Þ

ð10Þ

where sgn(x) function returns the sign of the parameter x. Although, truncation of the singular values using soft thresholding function gives filtering results better than the results obtained with traditional hard thresholding function but this truncation procedure also has some critical issues. Firstly, this procedure treats each singular value equally, and as a result soft thresholding operator in (9) reduces each singular value with the same amount which in turn will induce the deviation when the filtered results are reconstructed by the inverse SVD. Secondly, this procedure ignores the prior knowledge we often have on the matrix singular values. For instance, the larger singular values are generally associated with the important signal features, and thus they should better shrink less to preserve the significant image features. Clearly, this procedure in conjunction with soft thresholding operation fails to encash such prior knowledge. To resolve the issues raised in above truncation procedure, Gu et al. [18] recently proposed a new truncation procedure which generalizes the former and greatly improves its flexibility. Instead of using fixed threshold for the truncation of singular

Multimed Tools Appl

values, they obtained the appropriate weights w1, w2, …, wN corresponding to each singular value by using the prior knowledge about singular values. Then these weights are applied as thresholds for truncating the noisy singular values βi’s using soft thresholding function as: 0

αi ¼ T s ðβ i ; wi Þ ¼ sgnðβ i Þmaxðjβ i j−wi ; 0Þ

ð11Þ

The truncated singular values α′i, i = 1, …. N are the diagonal elements of the diagonal matrix Σ Xb . Now using (7), we can obtain the filtered matrix Xb . The filtering approach used in the proposed method, elaborated in the next section, is also inspired by the above suggested idea of using separate thresholds for all singular values based on their significance.

3 Proposed method Performing the SVD-based filtering in wavelet domain is not a new concept, but it has been used occasionally. Hou [19] applied an edge-adaptive SVD filtering in wavelet domain for image denoising. In this filtering method, initially each detail subband are divided into blocks and then the presence of an edge structure is checked in each block based on the median value of its singular values. Thereafter, the thresholds λedge and λnonedge, respectively, are applied on the edge-present blocks and edge-absent blocks according to (9) to obtain the truncated singular values for these blocks. The thresholds are computed by using the noise variance b σ 2noise which in turn is estimated through a robust median estimator [12]. The above filtering method has some issues that need to be resolved. Firstly, each of the above mentioned thresholds is used for multiple blocks (or matrices) to get their low rank approximations, while in Section 2, we have highlighted the problems with the case where even single threshold is used for filtering a single matrix. Secondly, the decision about the applicability of these thresholds on a particular block takes the local features of that block into consideration, but computation of any of these thresholds does not consider the local characteristics of the blocks. Thirdly, the term noise variance b σ 2noise used in the computation of thresholds is kept fixed through all resolution scales, while the noise strength decreases with the increment of resolution scale. Therefore, the noise variance is better to be estimated separately at each resolution scale. To resolve above mentioned issues, a new adaptive patch-based weighted-SVD filtering technique is proposed to suppress the noise from wavelet coefficients while preserving edges. This scheme includes the computation of noise level and estimation of weights (thresholds) for patch-based SVD filtering of wavelet coefficients. The main stages of the proposed denoising method are illustrated in Fig. 1. Suppose that a given spatial domain image f = {f(x, y), x = 1, …, M, y = 1, … N} has been corrupted by additive noise as: g ¼ f þn

ð12Þ

where g and f are the noisy image and the original image, respectively, n represents the noise and is modeled as i.i.d. Gaussian with zero mean and standard deviation σ, denoted as N ð0; σ2 Þ. The work presented in this paper deals with gray-scale images,

Multimed Tools Appl

noisy image

discrete wavelet transform

G

singular values truncation

inverse singular value decomposition

local noise variance estimation

non-local similar patches estimation

weight (threshold) vector estimation

singular value decomposition

aggregation

F

inverse wavelet transform

denoised image

Fig. 1 Proposed denoising method

which have a Gaussian distribution and hence there is no consideration of distribution of noise. Initially, the discrete wavelet transform W is applied on input noisy image g to obtain its decomposition G into different frequency subbands labeled as LLJ, LHj, HLj, and HHj, j = 1,2,….., J, where j is the scale and J is the coarsest scale. The larger j is, the coarser the scale is. G ¼ W ðg Þ

ð13Þ

Let us use the notations S1j , S2j and S3j for the detail subbands LHj, HLj and HHj, respectively. Consider the detail subband Slj, j = 1,…..,J and l = 1, 2, 3. We will proceed as follows. The local noise variance b σ 2noise; j can be estimated as: σ noise; j ¼ b

medianðjui jÞ ; ui ∈S 3j 0:6745

ð14Þ

The expression on right-hand side of Eq. (13) is a robust median estimator generally used to estimate noise variance from highest frequency subband, i.e., the subband S13) [12]. Note that the additive noise model for image g in spatial domain in Eq. (11) is also applicable for the subband Slj in wavelet domain. Thus, we have: S lj ¼ wlj þ nlj ;

ð15Þ

where the noisy subband Slj is obtained after adding the noise nlj in its noiseless counterpart wlj. The objective of the proposed adaptive patch-based weighted SVD j filtering method is to obtain w b l , the estimate of wlj, from its noisy observation Slj. For a local patch yk in the subband Slj, we can search for its non-local similar patches across the subband (for simplicity, in a large enough local search window) by methods such as block matching [10]. By stacking those non-local similar patches into a matrix, denoted by Yk, we have Yk = Wk + Nk, where Wk and Nk are the patch matrices of original noiseless subband wjl and noise njl, respectively. The main intuition is to obtain an estimate Wb k of patch matrix Wk from Yk, and the low rank matrix approximation method

Multimed Tools Appl

(i.e., the concept of truncated SVD) can be used for this purpose. Similar to (6), we have the SVD of Yk as: Y k ¼ U Y k Σ Y k V tY k

ð16Þ

We used a weight vector w = [w1, w2, …, wN], where the weight wi assigned to αi(Wk), the i-th singular value of Wk, can be computed as [18]: pffiffiffiffi 2 c Nσ ^ noise; j wi ¼ ðαi ðW k Þ þ εÞ

ð17Þ

where c > 0 is a constant, N is the number similar patches in Yk, and ε = 10− 16 is to avoid dividing by zero. As can be seen from Eq. (16) that the weight wi assigned to αi(Wk) is inversely proportional to αi(Wk), reason being the larger singular values that represent the energy of the major components of Wk should better be shrink less, so small weights should be associated with them. However, there is still a problem that the singular values αi(Wk) are not available. With the assumption that the noise energy is evenly distributed over each subspace by the basis pair of U Y k and V Y k , the initial αi(Wk) can be estimated as: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi αbi ðW k Þ ¼

2

max β 2i ðY k Þ−N σ ^noise; j ; 0 ;

ð18Þ

biðW k Þ is used in Eq. where βi(Yk) is the i-th singular value of Yk. The initial estimate α (16) to compute the weight wi. Note that the obtained weights w1, w2, …, wN are biðW k Þ are always sorted in a nondefinitely to be in a non-descending order since α ascending order. Instead of using the fixed value of threshold, the weight wi is used as threshold to obtain the truncated singular value α'i corresponding to βi(Yk). Therefore, Eq. (10) can be reused as: 0

αi ¼ T s ðβ i ; wi Þ ¼ sgnðβ i Þmaxðjβ i j−wi ; 0Þ

ð19Þ

The truncated singular values α'i, i = 1, …. N are the diagonal elements of the diagonal matrix Σ Wbk . SVD provides an explicit representation of the range and null space of a matrix X. The left singular vectors corresponding to the non-zero singular values of X span the range of X, while the right singular vectors corresponding to zero singular values of X span the null space of X. As the proposed approach applies the SVD on above defined matrix Yk. Although, initially Yk is a full-rank matrix, but after performing the proposed weighted SVD filtering we get a low-rank (say, r-rank) approximation of Yk. That is, the diagonal matrix Σ Wbk will have r non-zero singular values on its diagonal. Thus, for construction of noise-free results, we will use the matrix U Wbk which contains r left singular vectors from U Y k (i.e., vectors those span the range of Yk) and the matrix V Wbk

which contains the N − r right singular vectors from V Y k .

Multimed Tools Appl

Now similar to (7), Wbk can be estimated using the inverse singular value decomposition as: bk ¼ U Σ V t W W bk Wbk Wb

ð20Þ k

By applying the above procedures to each local patch yk and aggregating all patches j together, the estimate w b l can be obtained. The above procedure will be used to obtain the noiseless estimate for all the detail subbands at each level and the final thresholded result of complete image in wavelet domain can be obtained as: b¼ F

 b 1j ; w b 2j ; w b 3j ; j ¼ J ; J −1; …; 1g LL J ; w

ð21Þ

At last, the desired denoised image fb can be obtained by applying an inverse discrete wavelet transform W −1 on thresholded wavelet domain image Fb:

ð

−1 b FÞ fb ¼ W

ð22Þ

In practice, we can run several more iterations of all the stages to enhance the denoising outputs. The whole denoising algorithm is summarized in Table 1.

4 Experimental results To validate the proposed denoising method, a number of experiments were conducted on several real test images corrupted by simulated additive white Gaussian noise at

Table 1 Algorithm for proposed denoising method

Multimed Tools Appl

six different power levels σ [10, 14, 20, 25, 26, 31]. The proposed algorithm has been iterated over ten different noise realizations for each standard deviation and the results are averaged over these ten runs. The test set comprises images from standard grayscale image dataset (http://decsai.ugr.es/cvg/CG/base.htm), as well as well-known images such as Lena, cameraman and peppers. A subset of grayscale test images shown in Fig. 2 is considered in the subsequent discussions for testing the performance.

4.1 Performance metrics The standard simplified denoising problem is to find a reasonably good estimate fb of its unknown noise-free version f from the corresponding noisy version g. One of the

Fig. 2 Images used in the comparisons

(a)

Lena

(b) cameraman

(c) Barbara

(d)

man

(e)

(f)

peppers

boat

Multimed Tools Appl

common and simplest performance metric is the mean square error (MSE), defined as i2 2 1 1 XM XN h ^f ðx; yÞ f ð x; y Þ− ð23Þ ‖ f − ^f ‖ ¼ MSE ¼ x¼1 y¼1 M N M N where ‖ ⋅ ‖2 is the Euclidean norm. Most denoising applications use the peak signal-to-noise ratio (PSNR) as a performance metric, which is defined in decibels (dB) for 8-bit grayscale images as   2552 PSNR ¼ 10log10 ð24Þ MSE A higher value of PSNR would normally reflect the better performance of denoising algorithm. The MSE and PSNR are appealing performance metrics because they are easy to calculate and have clear physical meanings which are also mathematically convenient in the context of optimization. But they are not very well matched to perceive visual quality, i.e., they do not measure the visual features of resulting image directly. Apart from this, the images with large differences in their psychovisual quality can have similar MSE (or PSNR) values. To resolve these issues, the structural similarity index (SSIM) [44] was proposed as another metric to compare images which correlates more appropriately with the human perception. It maps two images into an index in the interval [−1, 1], where higher values are given to more similar pairs of images X and Y, calculated as SSIMðX ; Y Þ ¼

ð2μX μY þ C 1 Þð2σX Y þ C 2 Þ ðμ2X þ μ2Y þ C 1 Þðσ2X þ σ2Y þ C 2 Þ

ð25Þ

where μX, μY, σ2X and σ2Y are the averages and variances of X and Y, σXY is the covariance between X and Y, both C1 and C2 are predefined constants. To measure performance of an edge preserving denoising filter, Pratt’s figure of merit (FOM) [33] is a widely adapted metric. The value of this metric reflects that how much the detected edges in the filtered image coincide with ideal edges (that is, edges in the original image). In other words, the value of FOM gives an idea about how effectively edges of original image are preserved in the filtered image. This metric is defined as: FOM ¼

X ND 1 1 i¼1 maxðN I ; N D Þ 1 þ αd 2i

ð26Þ

where NI and ND are the numbers of ideal and detected edge pixels, respectively, di is the distance between an edge point and the nearest ideal edge pixel, α is an empirical calibration constant (often 1/9) used to penalize displaced edges. The value of FOM is a number in the interval [0, 1], where 1 represents the optimal value, that is, the detected edges coincide with the ideal edges.

4.2 Parameters setting We estimate a set of parameters used by the proposed method in section 3: wavelet transform and its number of resolution scales (decomposition levels) J, δ, c, iteration

Multimed Tools Appl

number K, and patch size. To train the parameters, a set of about 20 images, different from those shown in the comparisons, is used. The images for this set are chosen from a test database which contains 52 test images (49 images from the standard grayscale image dataset (http://decsai.ugr.es/cvg/CG/base.htm) and three well-known images Lena, cameraman, and peppers, respectively). Obviously, the performance of the proposed approach sensitive to the choices of parameter setting. In fact, the final value of all the parameters were set by comparing the results obtained through variation in parameters. Once the parameters are set, they are kept fixed throughout the comparisons to other methods. A set of stationary wavelets [29] from Symlet, Coeflet, Daubechies and Biothogonal families is tested for effectiveness. According to our experiments, Symlet-8 (sym8) provided better results than other wavelet bases. In addition, four resolution levels (i.e., J = 4) achieved the best results. Thus, for all other wavelet-based methods (SURELET, Besov, Bayes and WASVD), the Symlet-8 with four decomposition levels are used for fair comparison. The iterative regularization parameter δ and the parameter c are fixed to 0.1 and 2.8, respectively for all noise levels. The values of these parameters are set on the basis of a number of experiments performed on a set of images. The iteration number K and patch size are set based on noise level. For higher noise level, it is desirable to choose bigger patches and run more times the iteration. By experience, we set the value of K to 2, 3, and 5 for σ ≤ 20, 20 < σ ≤ 30, 30 < σ, respectively. Patch size is set to 5 × 5, 7 × 7, and 9 × 9, respectively on these noise levels.

4.3 Comparisons To assess the denoising effectiveness, the proposed method is compared to state-of-the-art methods, namely, SURELET [2], Besov-Ball projections (in short, Besov) [8] and Bayes [6], WASVD [19]. Tables 2, 3 and 4, respectively contain the PSNR (in dB), SSIM and FOM values of the denoised images relative to their original images using these methods. The best values amongst all the methods are highlighted with bold face. The results shown in tables demonstrate that the proposed method is superior to the methods Besov, Bayes and WASVD regarding the PSNR and SSIM measures but sometimes these methods are competing in regarding FOM measure. When we compared with the SURELET method, the proposed method is at times contested in respect to all available performance measures. In Fig. 3 which compares the PSNR performances of our method with other well-known denoising methods, starting from (a) to (f) it can be observed that the proposed method performs better in comparison to others in most of the cases. Though sometimes the LLSURE method (Fig. 3f) and SURELET method (Fig. 3d and f) performs marginally better then the proposed one. Figures 4, 5, 6 and 7 demonstrate the visual comparison among all the methods, including the proposed one in respect to the sample test images: Lena, Barbara, Cameraman and Peppers, respectively. The Bayes method tends to produce smoothened results in homogeneous regions. Nevertheless, certain features such as edges are affected. As the proposed denoising method employs an adaptive weighted SVD-based filtering scheme to threshold the wavelet coefficients, it is possible to observe that such adaptive thresholding, in conjunction with the noise level computation, effectively reduces noise while preserving features of the image. This effect can be better seen in Figs. 4(h) and 6(h). The SURELET method produces a similar result on edges. However, it can be perceived from Figs. 4, 5, 6 and 7, that the proposed method outperforms SURELET in homogeneous regions, producing smoother results. That can be clearly observed in the various smooth regions in Fig. 7.

Multimed Tools Appl Table 2 Performance of various methods as measured by PSNR Methods

Lena σ = 10 σ = 15 σ = 20 σ = 25 σ = 30 σ = 35 Cameraman σ = 10 σ = 15 σ = 20 σ = 25 σ = 30 σ = 35 Barbara σ = 10 σ = 15 σ = 20 σ = 25 σ = 30 σ = 35 Man σ = 10 σ = 15 σ = 20 σ = 25 σ = 30 σ = 35 Boat σ = 10 σ = 15 σ = 20 σ = 25 σ = 30 σ = 35 Peppers σ = 10 σ = 15 σ = 20 σ = 25 σ = 30 σ = 35

LLSURE

SURELET

Besov

Bayes

WASVD

Proposed

Method

Method

Method

Method

Method

Method

33.78 31.93 30.84 29.93 29.20 28.41

34.76 32.56 31.30 30.33 29.15 28.36

33.57 31.48 30.15 29.04 28.33 27.74

33.54 31.32 30.10 29.43 28.56 28.19

33.45 31.02 29.09 27.98 26.96 26.03

34.89 33.05 31.79 30.85 29.71 28.82

34.61 32.55 31.27 30.34 29.86 29.47

33.77 31.44 29.79 28.55 27.28 26.78

34.13 32.66 31.17 30.29 29.40 28.84

35.16 32.93 31.58 30.34 29.54 29.03

34.74 32.11 30.20 28.79 27.68 26.84

36.78 34.39 32.70 31.37 30.71 30.18

31.64 29.43 27.93 26.95 25.85 24.91

32.28 30.02 28.47 27.30 25.94 24.85

30.62 27.98 26.53 25.35 24.58 24.20

31.10 29.54 27.50 26.67 25.47 24.97

32.00 29.57 27.83 26.68 25.74 25.03

33.43 31.33 29.92 28.67 27.76 26.99

32.25 30.35 29.11 28.31 27.44 26.89

32.41 30.47 29.23 28.52 27.32 26.64

31.99 29.62 28.32 27.34 26.68 25.98

31.44 29.76 28.50 27.36 27.04 25.94

31.35 29.09 27.49 26.39 25.52 24.50

32.19 30.20 28.91 27.82 27.44 26.98

32.61 30.60 29.35 28.34 27.24 26.70

33.21 31.17 29.85 28.93 28.06 27.40

32.49 30.40 28.85 27.83 26.84 26.12

32.07 30.14 28.80 27.48 26.77 26.25

32.05 29.76 28.26 27.18 26.14 25.34

33.57 31.51 30.61 29.41 28.66 27.90

33.20 31.78 30.73 29.88 29.05 27.90

33.38 31.57 31.05 29.65 28.70 28.43

32.24 30.60 29.52 28.39 27.78 27.11

32.38 30.51 29.59 28.88 27.95 26.15

31.89 30.04 28.68 27.61 26.76 25.98

33.64 31.67 30.83 29.54 29.19 28.77

Multimed Tools Appl Table 3 Performance of various methods as measured by SSIM Methods LLSURE

SURELET

Besov

Bayes

WASVD

Proposed

Method

Method

Method

Method

Method

Method

Lena σ = 10

0.89

0.90

0.85

0.87

0.85

0.92

σ = 15

0.86

0.86

0.81

0.83

0.78

0.88

σ = 20

0.84

0.84

0.78

0.80

0.70

0.84

σ = 25

0.81

0.82

0.76

0.78

0.63

0.82

σ = 30

0.78

0.80

0.73

0.76

0.57

0.81

σ = 35

0.76

0.78

0.71

0.74

0.52

0.79

Cameraman σ = 10

0.92

0.93

0.83

0.89

0.87

0.94

σ = 15

0.89

0.90

0.81

0.84

0.78

0.90

σ = 20

0.80

0.87

0.78

0.80

0.69

0.88

σ = 25

0.84

0.85

0.77

0.76

0.62

0.86

σ = 30

0.81

0.83

0.75

0.74

0.55

0.83

σ = 35

0.77

0.81

0.73

0.59

0.49

0.81

Barbara σ = 10

0.91

0.90

0.87

0.85

0.88

0.92

σ = 15

0.86

0.85

0.80

0.79

0.81

0.88

σ = 20

0.81

0.80

0.76

0.74

0.74

0.85

σ = 25

0.80

0.75

0.71

0.70

0.68

0.80

σ = 30

0.74

0.71

0.67

0.67

0.62

0.78

σ = 35

0.70

0.67

0.65

0.64

0.57

0.74

Man σ = 10

0.88

0.88

0.87

0.86

0.83

0.88

σ = 15

0.83

0.83

0.81

0.80

0.75

0.85

σ = 20

0.79

0.79

0.77

0.75

0.68

0.79

σ = 25

0.75

0.75

0.73

0.71

0.61

0.77

σ = 30

0.73

0.73

0.70

0.69

0.56

0.73

σ = 35

0.70

0.70

0.67

0.67

0.51

0.71

Boat σ = 10

0.90

0.90

0.86

0.86

0.85

0.92

σ = 15

0.86

0.86

0.82

0.80

0.77

0.89

σ = 20

0.84

0.82

0.78

0.76

0.70

0.84

σ = 25

0.79

0.79

0.75

0.73

0.63

0.81

σ = 30

0.75

0.76

0.73

0.71

0.57

0.78

σ = 35

0.72

0.74

0.71

0.68

0.52

0.75

Peppers σ = 10

0.88

0.88

0.82

0.84

0.81

0.88

σ = 15

0.85

0.84

0.77

0.79

0.74

0.84

σ = 20

0.82

0.81

0.73

0.76

0.67

0.82

σ = 25

0.80

0.79

0.69

0.73

0.60

0.80

σ = 30

0.78

0.77

0.67

0.71

0.55

0.77

σ = 35

0.76

0.75

0.65

0.56

0.50

0.76

Multimed Tools Appl Table 4 Performance of various methods as measured by FOM LLSURE

SURELET

Besov

Bayes

WASVD

Proposed

Method

Method

Method

Method

Method

Method

Lena σ = 10

0.82

0.89

0.84

0.89

0.89

0.90

σ = 15 σ = 20

0.76 0.72

0.86 0.82

0.82 0.75

0.81 0.84

0.81 0.79

0.87 0.85

σ = 25

0.73

0.79

0.70

0.76

0.71

0.81

σ = 30

0.72

0.76

0.64

0.66

0.66

0.79

σ = 35

0.67

0.72

0.63

0.67

0.62

0.74

Cameraman σ = 10

0.80

0.90

0.92

0.91

0.91

0.92

σ = 15

0.75

0.84

0.85

0.87

0.88

0.88

σ = 20 σ = 25

0.72 0.65

0.80 0.76

0.68 0.64

0.74 0.62

0.74 0.64

0.86 0.79

σ = 30

0.67

0.74

0.67

0.59

0.65

0.74

σ = 35

0.61

0.70

0.60

0.57

0.61

0.71

Barbara σ = 10

0.89

0.92

0.90

0.91

0.92

0.93

σ = 15

0.79

0.87

0.84

0.87

0.88

0.92

σ = 20

0.77

0.85

0.80

0.83

0.80

0.91

σ = 25 σ = 30

0.75 0.71

0.78 0.76

0.73 0.70

0.78 0.72

0.76 0.71

0.87 0.86

σ = 35

0.66

0.73

0.69

0.69

0.69

0.82

Man σ = 10

0.79

0.88

0.90

0.91

0.92

0.90

σ = 15

0.76

0.84

0.86

0.88

0.87

0.87

σ = 20

0.72

0.80

0.81

0.84

0.81

0.82

σ = 25

0.71

0.79

0.79

0.80

0.75

0.80

σ = 30 σ = 35

0.69 0.65

0.78 0.72

0.75 0.70

0.76 0.74

0.74 0.70

0.79 0.74

Boat σ = 10

0.80

0.90

0.93

0.92

0.94

0.93

σ = 15

0.77

0.85

0.88

0.89

0.92

0.91

σ = 20

0.71

0.84

0.80

0.86

0.88

0.89

σ = 25

0.68

0.79

0.77

0.83

0.81

0.86

σ = 30

0.67

0.73

0.75

0.78

0.74

0.81

σ = 35 Peppers

0.66

0.75

0.69

0.75

0.70

0.80

σ = 10

0.84

0.88

0.82

0.85

0.87

0.88

σ = 15

0.81

0.85

0.73

0.81

0.78

0.84

σ = 20

0.78

0.81

0.70

0.80

0.69

0.82

σ = 25

0.76

0.80

0.61

0.75

0.69

0.81

σ = 30

0.74

0.78

0.63

0.65

0.60

0.78

σ = 35

0.72

0.72

0.56

0.52

0.56

0.73

Multimed Tools Appl

Fig. 3 PSNR performance graphs for test images

Multimed Tools Appl

(a) Original

(e) Besov

(b) noisy (σ =30)

(c) LLSURE

(f) Bayes

(g) WASVD

(d) SURELET

(h) Proposed

Fig. 4 Denoising results for image Lena

The Besov and Bayes methods fail to smoothen up images when noise increases to higher levels. These produce good results at lower σ values, but give poor denoised images at higher noise levels. The proposed method overcomes such problem as some crucial parameters are set based on noise levels.

(a) Original

(e) Besov

(b) noisy (σ =30)

(c) LLSURE

(d) SURELET

(f) Bayes

(g) WASVD

(h) Proposed

Fig. 5 Denoising results for image Barbara

Multimed Tools Appl

(a) Original

(e) Besov

(b) noisy (σ =30)

(c) LLSURE

(d) SURELET

(f) Bayes

(g) WASVD

(h) Proposed

Fig. 6 Denoising results for image Cameraman

The outcome of the proposed algorithm, along with the noise level, also depends on the local structures present in the given image. The ‘Peppers’ image is a special case; as this image has heterogeneous composition of elements in foreground and the background is totally missing as the foreground elements are dominating. Even though the proposed algorithm was able to get the smooth edges but it was not able to capture the depth of the groove which is

(a) Original

(b) noisy (σ =30)

(c) LLSURE

(d) SURELET

(e) Besov

(f) Bayes

(g) WASVD

(h) Proposed

Fig. 7 Denoising results for image Peppers

Multimed Tools Appl

having the shadow effect and the intensity values ranged in the black to gray range. This is still a challenge for the current formulation and also it can be seen that none of the standard algorithms used for comparative study were able to capture the same.

5 Conclusions This paper presented a new edge-preserving image denoising method in wavelet domain which employs an adaptive patch-based weighted-SVD filtering as a thresholding scheme. In this filtering method, we considered a group of non-local similar patches for each local patch in the subband under consideration and the weights are estimated for this patch group based on the significance of its singular values. These weights are then used as a threshold for SVD filtering of the patch group considered. Such procedure is applied to each local patch and all patches are aggregated together to obtain the thresholded results in the subband. The proposed method has several desirable features. First, the approach of obtaining the filtered outcome of a local patch by considering the effects of its non-local similar patches enhances the denoising performance. Second, the method shows its great ability to preserve fine details and geometrical structures in the original image, as it shrinks each singular value by an amount (threshold) which is estimated on the basis of the noise strength of that singular value. Third, estimating the term noise variance, used in the computation of weights (thresholds), locally at each resolution scale makes it more beneficial as it takes the noise strength at that scale into consideration. Both the quantitative and qualitative analysis of the results indicates that the proposed method effectively suppresses Gaussian noise without eliminating the important image details. Experiments demonstrate that the proposed method produces better results when compared with results obtained from other state-of-the-art denoising methods.

References 1. Antoniadis A, Fan J (2001) Regularization of wavelet approximations. J Am Stat Assoc 96(455):939–967 2. Blu T, Luisier F (2007) The SURE-LET approach to image denoising. IEEE Trans Image Process 16(11): 2778–2786 3. Buades A, Coll B, Morel JM (2005) A non-local algorithm for image denoising. In Proceeding of IEEE International Conference on Computer Vision and Pattern Recognition, vol. 2. San Diego, CA, USA: IEEE Press, pp. 60–65 4. Cai J-F, Cand’es EJ, Shen Z (2010) A singular value thresholding algorithm for matrix completion. SIAM J Optim 20(4):1956–1982 5. Chang S, Yu B, Vetterli M (2000) Spatially adaptive wavelet thresholding based on context modeling for image denoising. IEEE Trans Image Process 9(9):1522–1531 6. Chang S, Yu B, Vetterli M (2000) Adaptive wavelet thresholding for image denoising and compression. IEEE Trans Image Process 9(9):1532–1546 7. Chipman H, Kolaczyk E, McCulloc R (1997) Adaptive Bayesian wavelet shrinkage. J Am Stat Assoc 440(92):1413–1421 8. Choi H, Baraniuk RG (2004) Multiple wavelet basis image denoising using Besov ball projections. IEEE Signal Process Lett 11(9):717–720 9. Dabov K, Foi A, Katkovnik V, Egiazarian K (2006) Image denoising with block-matching and 3D filtering. SPIE Electron Imaging Algorithm Syst 6064:606414–1–606414–12 10. Dabov K, Foi A, Katkovnik V, Egiazarian K (2007) Image denoising by sparse 3-d transform-domain collaborative filtering. IEEE Trans Image Process 16(8):2080–2095 11. Donoho DL (1995) De-noising by soft-thresholding. IEEE Trans Inf Theory 41(3):613–627 12. Donoho DL, Johnstone IM (1994) Ideal spatial adaptation via wavelet shrinkage. Biometrika 81:425–455

Multimed Tools Appl 13. Donoho DL, Johnstone IM (1995) Adapting to unknown smoothness via wavelet shrinkage. J Am Stat Assoc 90(432):1200–1224 14. Gao HY (1998) Wavelet shrinkage denoising using the nonnegative garrote. J Comput Graph Stat 7(4):469–488 15. Gao HY, Bruce AG (1997) WaveShrink with firm shrinkage. Stat Sin 7:855–874 16. Golub GH, Van Loan CF (1983) Matrix computations. John Hopkins University, Press, Baltimore 17. Gonzalez RC, Woods RE (2008) Digital image processing. Prentice-Hall, Upper Saddle River 18. Gu S, Zhang L, Zuo W, Feng X (2014) Weighted nuclear norm minimization with application to image denoising. IEEE Conf Comput Vis Pattern Recogn 19. Hou Z (2003) Adaptive singular value decomposition in wavelet domain for image denoising. Pattern Recogn J Pattern Recogn Soc 36:1747–1763 20. Jain P, Tyagi V (2013) Spatial and frequency domain filters for restoration of noisy images. IETE J Educ 54(2):108–116 21. Jain P, Tyagi V (2014) A survey of edge-preserving image denoising methods. Inf Syst Front 1–12. doi: 10. 1007/s10796-014-9527-0 22. Jain P, Tyagi V (2015) An adaptive edge-preserving image denoising technique using tetrolet transforms. Vis Comput 31(5):657–674. doi:10.1007/s00371-014-0993-7 23. Jain P, Tyagi V (2015) LAPB: locally adaptive patch-based wavelet domain edge-preserving image denoising. Inform Sci 294:164–181. doi:10.1016/j.ins.2014.09.060 24. Konstantinides K, Natarajan B, Yovanof GS (1997) Noise estimation and filtering using block-based singular value decomposition. IEEE Trans Image Process 6(3):479–483 25. Konstantinides K, Yao K (1988) Statistical analysis of effective singular values in matrix rank determination. IEEE Trans Acoust Speech Signal Process 757–763 26. Luo G (2006) Fast wavelet image denoising based on local variance and edge analysis. Int J Intell Technol 1(2):165–175 27. Maggioni M, Katkovnik V, Egiazarian K, Foi A (2013) Nonlocal transform-domain filter for volumetric data denoising and reconstruction. IEEE Trans Image Process 22(1):119–133 28. Nason GP (1996) Wavelet shrinkage by cross-validation. J R Stat Soc B 58:463–479 29. Nason GP, Silverman BW (1995) The stationary wavelet transform and some statistical applications. In: Lecture notes in statistics: wavelets and statistics, Springer-Verlag, Berlin, pp. 281–300 30. Orchard J, Ebrahimi M, Wong A (2008) Efficient nonlocal-means denoising using the SVD. In: IEEE international conference on image processing. San Diego, CA, USA 1732–1735 31. Perona P, Malik J (1990) Scale-space and edge detection using anisotropic diffusion. IEEE Trans Pattern Anal Mach Intell 12(7):629–639 32. Portilla J, Strela V, Wainwright M, Simoncelli E (2003) Image denoising using scale mixtures of Gaussians in the wavelet domain. IEEE Trans Image Process 12(11):1338–1351 33. Pratt WK (2012) Digital image processing. Wiley, New York 34. Qiu T, Wang A, Yu N, Song A (2013) LLSURE: local linear sure-based edge-preserving image filtering. IEEE Trans Image Process 22(1):80–90 35. Rudin L, Osher S D1994] Total variation based image restoration with free local constraints. Proc IEEE Int Conf Image Process 1:31–35, Austin, Texas, USA 36. Rudin L, Osher S, Fatemi E (1992) Nonlinear total variation based noise removal algorithms. Phys D 60:259–268 37. Sendur L, Selesnick IW (2002) Bivariate shrinkage with local variance estimation. IEEE Sig Process Lett 9(12):438–441 38. Shao L, Yan R, Li X, Liu Y. From heuristic optimization to dictionary learning: a review and comprehensive comparison of image denoising algorithms 1–14 39. Shapiro L, Stockman G (2001) Comput Vis. Prentice Hall 40. Silva RD, Minetto R, Schwartz WR, Pedrini H D2012] Adaptive edge-preserving image denoising using wavelet transforms. Pattern Anal Appl. doi:10.1007/s10044-012-0266-x, Springer-Verlag 41. Tomasi C, Manduchi R (1998) Bilateral filtering for gray and color images. In proceeding of 6th International Conference on Computer Vision., Bombay, India 839–846 42. Van Der Veen AJ, Deprettere EF, Lee Swindlehurst A (1993) Subspace-based signal analysis using singular value decomposition. Proc IEEE 81:1277–1308 43. Vidakovic B (1998) Nonlinear wavelet shrinkage with Bayes rules and Bayes factors. J Am Stat Assoc 93(441):173–179 44. Wang Z, Bovik A, Sheikh H, Simoncelli E (2004) Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process 13(4):600–612 45. Weyrich N, Warhola GT (1998) Wavelet shrinkage and generalized cross validation for image denoising. IEEE Trans Image Process 7(1):82–90 46. Wongsawat Y, Rao K, Oraintara S (2005) Multichannel SVD based image denoising. IEEE Int Symp Circuits syst 6:5990–5993

Multimed Tools Appl

Paras Jain received his M.E. degree in Computer Engineering from SGSITS, Indore, India, affiliated with Rajiv Gandhi Technical University Bhopal, India in 2010. He has obtained his Ph.D. degree in CSE under the guidance of Dr. Vipin Tyagi from Jaypee University of Engineering and Technology, Guna, INDIA.

Dr. Vipin Tyagi, FIETE is faculty in CSE at Jaypee University of Engg and Technology, Raghogarh, Guna (MP) India. He has about 20 years of teaching and research experience. He was President of Engineering Sciences Section of the Indian Science Congress Association for the term 2010–11, and recorder for the term 2008–2010. He is a Life Fellow of the Institution of Electronics and Telecommunication Engineers. He is a senior life member of Computer Society of India. He is actively associated with such professional societies. He is member of Academic-Research and Consultancy committee of Computer Society of India. He is head of CSI Special Interest Group on Cyber Forensics. He is life member of CSI, Indian Remote Sensing Society, CSTA, ISCA IEEE, and International Association of Engineers. He was nominated to visit Czech Republic from INSA, New Delhi for 2 weeks in May 2012. He has published more than 50 papers in various journals, advanced research series and has attended several national and international conferences in India and abroad. He is Principal Investigator of research projects funded by DRDO, MP Council of Science and Technology and CSI.