A design framework for hybrid approaches of image noise estimation and its application to noise reduction

J. Vis. Commun. Image R. 23 (2012) 812–826

Contents lists available at SciVerse ScienceDirect

J. Vis. Commun. Image R. journal homepage: www.elsevier.com/locate/jvci

A design framework for hybrid approaches of image noise estimation and its application to noise reduction Shih-Ming Yang ⇑, Shen-Chuan Tai Department of Electrical Engineering, Institute of Computer and Communication Engineering, National Cheng Kung University, No. 1, University Rd., Tainan 701, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 9 June 2011 Accepted 16 April 2012 Available online 27 April 2012 Keywords: Noise estimation Noise reduction Image restoration Image enhancement Bilateral filter Non-local mean filter Laplacian operator Sobel operator

a b s t r a c t Noise estimation is an important process in digital imaging systems. Many noise reduction algorithms require their parameters to be adjusted based on the noise level. Filter-based approaches of image noise estimation usually were more efficient but had difficulty on separating noise from images. Block-based approaches could provide more accurate results but usually required higher computation complexity. In this work, a design framework for combining the strengths of filter-based and block-based approaches is presented. Different homogeneity analyzers for identifying the homogeneous blocks are discussed and their performances are compared. Then, two well-known filters, the bilateral and the non-local mean, are reviewed and their parameter settings are investigated. A new bilateral filter with edge enhancement is proposed. A modified non-local mean filter with much less complexity is also present. Compared to the original non-local mean filter, the complexity is dramatically reduced by 75% and yet the image quality is maintained. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Noise can be introduced into digital images during the acquisition and transmission processes. Estimation of the noise level is very important for a variety of digital image and video processing algorithms, such as denoising [1–8], super-resolution [9], feature extraction [10], and motion estimation [11]. The estimated results are usually used to adapt the algorithm parameters to the amount of noise instead of fixed values. The noise model which was often assumed is the additive white Gaussian noise (AWGN) process given by:

Iðx; yÞ ¼ f ðx; yÞ þ gðx; yÞ;

ð1Þ

where x and y are the vertical and horizontal coordinates of a pixel, I(x, y), f(x, y) and g(x, y) are the observed noisy image, the original ideal image and the additive Gaussian noise respectively. The image has width W and height H, and each pixel has an intensity value between 0 and 255. The goal of image noise estimation is to estimate the standard deviation rn of the Gaussian noise while the only information available is the noisy image I(x, y). The deviation rn is referred to as the noise level. Many algorithms for estimating the AWGN deviation in images have been proposed. They can be classified into wavelet-based [13,14], fuzzy-based [15], filter-based [16–18], block-based ⇑ Corresponding author. Fax: +886 6 2388673. E-mail addresses: [email protected] (S.-M. Yang), [email protected] (S.-C. Tai). 1047-3203/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jvcir.2012.04.007

[12,19,20], or hybrid [21,22]. The challenge is to separate noise from image features or details. In filter-based approaches, the noisy image is first filtered by a low-pass filter to remove the noise. Then the noise variance is computed from the difference between the noisy image and the filtered image. The main difficulty of filter-based approach is that the difference image is assumed to be the noise but this assumption is not always true, especially for images with fine structures or details. In blocked-based approaches, the noisy image is tessellated into a number of blocks. The noise variance is then computed based on a set of homogeneous blocks. The main issue of block-based approaches is how to identify the homogeneous blocks. A new hybrid approach was proposed in [22]. The noisy image is first tessellated into blocks. To identify the homogeneous blocks and exclude structures or details from contributing to the estimation of noise variance, a Sobel edge detection operator with a self-determined threshold is applied to each block. Then a filter operation, followed by an averaging of the convolutions over the selected blocks, provides a very accurate estimation of noise variance. It has been shown that the new hybrid approach achieves better performance in terms of accuracy, reliability and stability compared to other recently published methods. The goal of image noise estimation is to provide the estimated noise variance for the image processes that follow. Besides the noise variance, some side information generated during the estimation process could be valuable, too. For example, the Sobel edge detection operator can provide the edge information. But, not all image processes require the same side information. Therefore, it is of

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

interest to know the possibility of using homogeneity analyzers other than the Sobel operator in the noise estimation process. The main contributions of this work are the follows. First, a design framework for hybrid approaches of image noise estimation is proposed. Then, the performances of four homogeneity analyzers in terms of estimation accuracy and reliability are evaluated. Finally and most importantly, information provided by noise estimation is successfully integrated with two well-known non-linear noise reduction filters, the bilateral filter and the non-local mean filter. A new bilateral filter equipped with edge enhancement is proposed. A modified non-local mean filter with less complexity is also present. All of these improvements are made possible by the side information provided by noise estimation process. The remaining of this paper is organized as follows: Section 2 presents the proposed design framework. Section 3 introduces four homogeneity analyzers and Section 4 compares their performances in terms of estimation accuracy and reliability. Section 5 reviews the bilateral filter and investigates its parameter settings. A modified bilateral filter with edge enhancement is also proposed there. Section 6 reviews the non-local mean filter and investigates its parameter settings. A modified version of the non-local mean filter with much less complexity is presented there, too. Section 7 concludes the paper. 2. A design framework for hybrid approaches of image noise estimation Based on our previous work [22], a design framework for combining the hybrid approach of image noise estimation and noise reduction is presented in Fig. 1. The intention is to identify the homogeneous blocks first. Then, it will become much easier to separate noise from image itself on these blocks. The entire system begins with a block-based operation. The noisy image is tessellated into adjacent non-overlapping blocks with the block size equal to Wb. The choice of Wb has a great impact on the performance of the noise estimator. Considering performance and complexity, the best choice of Wb was shown to be 5 [22]. Let Bxy denote the block centered at (x, y) (2).

Bxy ¼

     Wb Wb and jn  yj 6 ; ðm; nÞjjm  xj 6 2 2

bc is the floor function:

ð2Þ

The homogeneity measure Hxy of each block is obtained by applying convolution with a homogeneity analyzer FH which is

consisted of Nf masks F Hi . The homogeneity measure is then defined as:

Hxy ¼

Nf X jF Hi  Bxy j i¼1

F Hi

ð3Þ

X

 Bxy ¼

wðm; nÞIðm; nÞ;

ðm;nÞ2Bxy

where w(m, n)’s are the coefficients of the masks. Four different operators are proposed to be the homogeneity analyzer and will be discussed in the next section. Depending on the choice of the operator, Nf could be 1, 2, 4, or 8. The main challenge of block-based approaches is to decide a threshold value for the homogeneity measures and identify the homogeneous blocks. Since the only information available is the homogeneity measure Hxy of each block, a self-determined threshold can be obtained based on the statistics of these measures. The histogram of Hxy is computed first. Then, the threshold value Hth is selected to be the Hxy value when the accumulated histogram reaches P% of the whole image. An image block is claimed to be homogeneous if the homogeneity measure Hxy is less than the threshold value of Hth. Considering the complexity and the estimation performance, the best choice of P was shown to be 10 [22]. After identifying the homogeneous blocks, each block is filtered by the Laplacian filter (4), which is the difference between two second-order masks, to reduce the interference of noise and exclude the image structures from contributing to the estimation [16].

3 2 3 2 3 1 2 1 0 1 0 1 0 1 7 6 7 6 6 7 2 5 LA ¼ 4 0 4 0 5  2  4 1 4 1 5 ¼ 4 2 4 1 2 1 0 1 0 1 0 1 2

rest ¼ BH

p1=2 2

1 6N H

X jLA  Bxy j BH

Finally, the estimated noise variance and the side information generated by the homogeneity analyzer are integrated into the noise reduction process. There are only two parameters in our estimation algorithm, the block width (Wb) and the percentage (P). P controls how many

Identifying Homogeneous Blocks

convolution

Side information convolution

Restored Image

Removing Image Noise

Noise Variance

ð5Þ

¼ fBxy jHxy < Hth g

Finding Homogeneity Measures

Tessellating Image

ð4Þ

Let BH denote the set of homogeneous blocks and BH denote the total number of pixels in BH . The estimated noise deviation is shown [16] to be given by (5)

Homogeneity Analyzer Noisy Image

813

Averaging Convolutions

Filtering Homogeneous Blocks

Fig. 1. Design framework for image noise estimation and reduction.

814

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

blocks involved in the final averaging process. It was reported in [22] that better performance could be achieved with smaller block size and lower percentage at low noise levels because the interferences by image features could be avoided. On the other hand, both parameters should be set to larger values at high noise levels because more samples will provide more accurate and reliable estimation. For a single image, since there was no prior information available, both parameters could only be set to fixed values. For video sequences, the estimation result of the previous frame could provide the valuable information for the next frame [23]. The proposed design framework is inspired by our previous work in [22] where the Sobel operator was applied as the homogeneity analyzer. The performance comparisons and complexity analysis of our estimation algorithm against other approaches were reported there. It is assumed that the image intensity is represented by an integer of range 0–255. Pixel gray levels may be driven into saturation by noise. The design issue regarding to the clipping or saturation effect was also addressed there. Three quality measures were defined to show the accuracy, reliability and stability of noise estimation algorithms. Our proposed algorithm was shown to achieve the best performance in terms of all these three measures and maintain low complexity. Here, we will investigate the performance the proposed framework using different homogeneity analyzers.

3. Homogeneity analyzers The key to the success of a block-based noise estimator is to identify the homogeneous area of the noisy image. Any operator which defines its own measure of homogeneity can be a candidate of homogeneity analyzers. Here, we consider four well-known operators which are Sobel, Laplacian, DSA (Directional Structure Analyzer), and LoG (Laplacian of Gaussian). 3.1. The Sobel operators (SBL) The first-order gradient is the most well-known measure for edge detection [24]. That’s why we chose the Sobel operators as the homogeneity analyzers in our previous work [22]. The Sobel masks with block size Wb = 5 are defined by (6) where Sv and Sh provide gradients in the vertical and horizontal directions, respectively.

-1

2

-1

4

-1

-1

axy ¼ tan1

-1

-1

4

-1

-1

-1

-1

Ld5

Ld6

ð7Þ

 Sv ; axy 2 ½p=2; p=2 Sh

ð8Þ

3.2. The Laplacian operators (LPN) The Laplacian operators generally are not used for edge detection because of two reasons [24]. First, they are sensitive to noise. Secondly, they are not able to detect the edge direction. But, for the purpose of homogeneity analysis, the second-order gradients provided by the Laplacian operators still can reflect the strength of local edge structure. Besides, the sign of the second derivative can tell us whether an edge pixel lies on the dark or light side of an edge. This information could be integrated with a noise reduction filter to improve the visual quality. The Laplacian operators with block size Wb = 5 are defined by (9). The homogeneity measure using the Laplacian operators is defined by (10).

2

0 6 60 6 LN ¼ 6 6 1 6 40 0

3

2 1 1 7 6 1 1 1 0 7 1 0 6 7 6 7 6 1 12 1 1 7 LD ¼ 6 0 0 7 6 1 1 1 0 5 4 1 0 0 1 0 0 1 1 0

1 0

0

0

1 1

3

7 1 7 7 12 0 0 7 7 7 0 0 1 5 0 1 1 0

0

ð9Þ HLPN xy ¼ jLN  Bxy j þ jLD  Bxy j

ð10Þ

-1

-1 -1

-1 4

4 -1

-1 -1

Ld2

-1

4 8

In addition to the gradients that detect the presence of edges, the direction of the gradient vector defined by (8) also tells us the direction of an edge. This is the information which can be integrated with noise filtering.

-1

-1

3 1 7 2 7 7 0 0 0 7 7 7 12 8 2 5 6 4 1 6 12

HSBL xy ¼ jSv  Bxy j þ jSh  Bxy j

-1

4

0 2 0 8

ð6Þ

4

Ld1

2 8

The homogeneity measure using the Sobel operators is defined by (7).

-1 -1

3 2 1 1 4 7 6 4 7 8 62 7 6 60 S 12 0 12 6 7 ¼ 0 h 7 6 7 6 8 0 8 4 5 4 2 8 1 2 0 2 1 1 4

1 6 64 6 Sv ¼ 6 66 6 44

-1

-1

-1

Ld3

Ld4

-1

-1

-1

-1

4

4

Ld7

Ld8

Fig. 2. Directional structure analyzers for Wb = 5.

-1

-1

815

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

3.3. The Directional Structure Analyzer (DSA) Amer and Dubois [19] proposed a structure-oriented noise estimator based on high-pass operators that measure homogeneity in eight different directions. The coefficients for each direction is {1, . . . , Wb  1, . . . 1}. An example for Wb = 5 is shown in Fig. 2. The homogeneity measure using DSA is defined by (11) where Ldi is illustrated in Fig. 2.

HDSA xy ¼

8 X jLdi  Bxy j

The digital approximation to the continuous LoG function is not unique. A 5  5 mask which can capture the essential shape of LoG is given by (12) [24,25]. The homogeneity measure using the LoG operator is given by (13).

2 LLoG

ð11Þ

0

0

1 0

0

3

7 6 1 2 1 0 7 60 7 6 7 ¼6 6 1 2 16 2 1 7 7 6 1 2 1 0 5 40 0 0 1 0 0

ð12Þ

i¼1

It has been shown [19] that DSA could provide accurate information about the image structure. In their algorithm, HDSA xy is only used as the preliminary criteria to identify three most homogeneous blocks and then decide the threshold value on local variances. The problem of using local variances as the measure for identifying the homogeneous area has been reported in [12,22]. Here, we will see how well DSA fits into our estimation algorithm. 3.4. Laplacian of Gaussian (LoG) As mentioned above, the Laplacian operator generally is not used in its original form for edge detection because it is sensitive to noise. The LoG operator, which combines the Laplacian second-order derivative and Gaussian smoothing, make it more robust to noise compared to other derivative operators [24].

HLoG xy ¼ jLLoG  Bxy j

ð13Þ

4. Performance comparisons of homogeneity analyzers 4.1. Experiment setup and quality measures For the Monte–Carlo simulations, 20 test images (Fig. 3) are used here. These test images are considered to be the original clean images. Computer generated Gaussian random noises are added to the clean images. The simulations were performed at noise levels (noise deviations) of range 1–30. The simulation for each image at each noise level is repeated 20 times. The quality measures for performance comparisons are the average estimation error En, the average variance of estimation errors r2En , and the Monte–Carlo variance V nM defined by (14).

Airfield

Mandrill

Barbara

Boat

Boots

Cameraman

Family

Fruits

Goldhill

Gray127

Harbour

Jet

Lena

LivingRoom

Lobster

Macbeth

Newspaper

Pepper

Toys

Wine

Fig. 3. Test images for performance comparisons of noise estimations.

816

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826 M X

Ei;n ¼ M1

i;n jri;n est ðmÞ  radd ðmÞj;

m¼1 NI X En ¼ N1I Ei;n ; i¼1

r2En ¼ N1I

NI X ðEi;n  En Þ2 :

ð14Þ

i¼1 1 ri;n est ¼ M

M X

ri;n est ðmÞ

m¼1

V nM ¼ N1I

NI M X X i;n 2 1 ðri;n est ðmÞ  rest Þ M i¼1

m¼1

2. NI is the total number of test images (Fig. 3). 3. Ei,n is the average estimation error of the M Monte–Carlo simulations for image i at noise level n. 4. En is the overall average estimation error at noise level n. It represents the accuracy of the noise estimator. For a good estimator, it should be as small as possible. 5. r2En is the overall average variance of estimation errors at noise level n. It represents the reliability of the noise estimator. For a good estimator, it should be as small as possible. 6. V nM is the variances of M repetitions averaged over all test images. It represents the stability of the noise estimator. It was first proposed in [22]. For a good estimator, it should be as small as possible.

Here, the subscript ‘‘est’’ denotes the estimated Gaussian noise deviation, the subscript ‘‘add’’ denotes the deviation of the simulated Gaussian noise added to the test images, the super- and subscripts i represent the index of the test images, and n denotes the simulated noise level. Also, 1. For each image at each noise level, the Monte–Carlo simulation is repeated M times and m is the index or counter of these repetitions. M was set to 20.

Fig. 4. Average estimation errors of DSA on single images using two, four, or eight directional masks.

Fig. 5. Overall Performances of DSA using two, four, or eight directional masks. (a) Average estimation errors. (b) Average variance of estimation errors. (c) Monte– Carlo variance.

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

4.2. Performance comparisons 4.2.1. Performances of DSA The DSA operators proposed in [19] required eight masks corresponding to eight directions. However, the other three homogeneity analyzers discussed here only take one or two masks. It is of interest to investigate the performances of DSA considering only two (Ld1 and Ld2) or four directions (Ld1–Ld4) (Fig. 2). The comparisons of DSA using two (DSA_D2), four (DSA_D4) and eight (DSA_D8) masks in terms of average estimation errors on single images are shown in Fig. 4. As it can be seen, DSA using eight masks only achieves better performances for images with fine details and for low noise deviations. The DSA masks are designed to catch image structures at different directions. There is no doubt that using more masks will perform better for images like Mandrill. But, for images like Lena, the advantage is not so obvious because there are not many fine details. When the noise deviations are high, the image structures are destroyed and the advantage of eight masks is reduced. The overall performances of DSA using two, four and eight masks are shown in Fig. 5. It can be concluded that DSA using two directional masks (DSA_D2) achieves lower estimation errors and better stability. We will compare DSA_D2 against the other homogeneity analyzers in the following.

Fig. 6. Performance comparisons on single images between different homogeneity analyzers for Wb = 5.

817

4.2.2. Comparisons of homogeneity analyzers The performance comparison in terms of average estimation errors on single image is shown in Fig. 6. ‘‘Lena’’ is the most popular test image for image processing algorithms. ‘‘Mandrill’’ is rich of details. ‘‘Barbara’’ has a lot of patterns. It can be observed that the Sobel operators provide the most accurate estimations for different types of images. The other operators perform equally well. The overall performance comparison in terms of average estimation errors or accuracy is shown in Fig. 7(a). The Sobel operators achieve the best overall performance. The Laplacian operators, DSA_D2 and the LoG operator perform equally well. Although the estimation errors of these three operators are higher at high noise deviations, the estimation errors are still lower than 5% of the noise deviations. If the parameters of the noise reduction filters are not too sensitive to the noise deviations, they still make good noise estimators.

Fig. 7. Comparisons of overall performances between different homogeneity analyzers for Wb = 5, (a) average estimation errors, (b) average variance of estimation errors, (c) Monte–Carlo variance.

818

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

The overall performance comparison in terms of average variance of estimation errors and Monte–Carlo variances are shown in Fig. 7(b) and (c). As it can be seen, all homogeneity analyzers perform equally well in terms of reliability and stability. Noted that the LoG operator has only one mask, it achieves the same performance with lower complexity compared to the others. In a short conclusion, we propose a design framework for hybrid approaches of image noise estimation. All the homogeneity analyzers discussed here fit into this framework and make accurate and reliable noise estimations. Especially, the best performance is achieved when the Sobel operator is adapted as the homogeneity analyzer. Its performance against other techniques was studied thoroughly and presented in [22]. In an image denoising application, the choice of homogeneity analyzers really depends on the side information which the noise reduction filter requires. 5. Bilateral filter 5.1. The filter formulation The bilateral filter [3] is a non-linear filter that smoothes images while preserving edges. Its application can be found in many image processing systems. A gentle and thorough introduction was presented in [26]. Let g^ðx; yÞ be the restored image, the output of the bilateral filter at pixel (x0, y0) can be formulated as follows.

g^BL ðx0 ; y0 Þ ¼

X 1 C BL ðx0 ; y0 Þ ðx;yÞ2Xðx

kBL ðx; y; x0 ; y0 ÞIðx; yÞ

8 > < W f 2 f3; . . . ; W max g rd 2 f1; . . . ; W f g > : rr ¼ r  radd ; r 2 f1; 2; 3; 4; 5g

; W max

8 rn 6 5 5 > > > <9 5 < rn 6 10 ¼ > 15 10 < rn 6 15 > > : 21 rn > 15 ð20Þ

The PSNR values between the original clean images and the restored images are recorded. The following results are observed. 5.2.2. Parameter settings for bilateral filter We begin with the selection of rd and rr. The performances of the bilateral filter using the test image Lena contaminated by AWGN with radd = 10 are shown in Fig. 8. It can be seen that: 1. Different selections of rd do not result in notable improvement or degradation when the other parameters are fixed. One half of the window size, Wh (18), is a good choice for rd. 2. Different selections of rr do make notable differences on PSNR values when the other parameters are fixed. In general, higher rr results in higher PSNR and blurred edges. 2rest or 3rest is a safe choice for rr, i.e., 2 or 3 for r (20). 3. As far as the window size (Wf) concerns, Wf = 3 is a good choice for low noise levels and Wf = 5 or 7 for medium to high noise levels.

ð15Þ

0 ;y0 Þ

where Xðx0 ; y0 Þ defined by (16) is the filtering window and kBL( x, y; x0, y0) defined by (17) is the weight that the pixel at (x,y) will contribute to the restoration of the pixel at (x0, y0). Wh, defined by (18), is one half of Wf which represents the window size of filtering operation.

Xðx0 ; y0 Þ ¼ fðx; yÞjjx  x0 j 6 W h and jy  y0 j 6 W h g kBL ðx; y; x0 ; y0 Þ ¼ exp 

ðx  x0 Þ2 þ ðy  y0 Þ2 2r2d

! exp 

ð16Þ

ðIðx; yÞ  Iðx0 ; y0 ÞÞ2 2r2r

!

(a)

ð17Þ W f ¼ 2W h þ 1

ð18Þ

CBL(x0, y0) defined by (19) is a normalization factor.

C BL ðx0 ; y0 Þ ¼

X

kBL ðx; y; x0 ; y0 Þ

ð19Þ

ðx;yÞ2Xðx0 ;y0 Þ

5.2. Searching for the best parameter settings 5.2.1. Experiment setup As it can be seen in (17), there are two parameters which control the behavior of the bilateral filter. They are rd and rr. In order to find out the relationship between the best parameter settings and the noise level (noise deviation) or image structure, we did exhaustive search. The same empirical study on the optimal parameter selection was done before [27–29]. However, in their work, the window size of the filter was not considered as one of the controlled parameters. They only considered the mean squared error (MSE) as the quality measure. Their conclusions suggested that rd is fixed to 2 and rr is set to nearly 2 times of the estimated noise deviation. We will present more detailed discussions and propose a better rule of thumb. Computer generated AWGN noise is added to standard test images like Lena, Mandrill, etc. (Fig. 3). The bilateral filter is applied to noisy image with the following parameter combinations.

(b)

(c) Fig. 8. Performance comparisons of the bilateral filter with different parameter settings.

819

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

4. A general rule of thumb will be: 8 ½W f ; rd ; rr  ¼ ½3; 1; rest  > < Low noise levelðrest < 5Þ Medium noise levelð5 6 rest 6 10Þ ½W f ; rd ; rr  ¼ ½3; 1; 3rest or½5; 2; 3rest  > : ½W f ; rd ; rr  ¼ ½5; 2; 3rest or½7; 3; 3rest  High noise levelðrest > 10Þ

ð21Þ The above discussions are based on PSNR values. If we consider noise reduction as the first priority, we can always choose larger Wf and rr. The price we pay will be more severe blurring. An example of the trade-off between smoothing noise and blurring is shown in Fig. 9. From Fig. 9(b) to (d), Wf is increased from 3 to 7. The image in Fig. 9(d) is restored with the largest Wf and achieves the best visual quality (by noise reduction), but it loses some details. On the other hand, the image in Fig. 9(b) is restored with the smallest Wf and preserves the edges as well as noise. The original intention of the bilateral filter is to preserve edges while smoothing images. Since it is still an averaging filter, blurring will always be inevitable. Therefore, we should always choose the parameter settings which achieve the best noise reduction results, and incorporate enhancement techniques to deal with the blurring. 5.2.3. Proposed modified bilateral filter with sharpening offset The way that the bilateral filter preserves edges is by using the range filter to enforce the gray level similarity. The range filter essentially is a Gaussian filter controlled by rr. It preserves edges by putting more weights on the neighbor pixels whose gray levels are closer to the central pixel. In order to enhance or sharpen edges, we need to do better than that. For example, if the central

pixel was on the light side of an edge, in order to enhance the edge, we should put more weights on the pixels of the light side and fewer weights on the pixels of the dark side. In other words, we need to shift the center of the range filter when we process edge pixels. It can be done by introducing an offset s(x0, y0) into the range filter (22), and the offset value should be adjusted based on some side information about the central pixel itself.

kBL ðx; y; x0 ; y0 Þ ¼ exp 

ðx  x0 Þ2 þ ðy  y0 Þ2 2r2d

 exp 

!

ðIðx; yÞ  Iðx0 ; y0 Þ  sðx0 ; y0 ÞÞ2 2r2r

! ð22Þ

The same idea of adding a sharpening offset was proposed in [7], too. They first use the LoG (Laplacian of Gaussian) operator as the pixel classifier. The offset and rr are adaptively selected based on the results of classification. The parameters involved in their algorithm eventually were decided by a training technique. We propose a new approach where all parameters are decided by heuristic. First of all, we choose the Laplacian operator (9) as our homogeneity analyzer for the noise estimation process. As mentioned in Section 3.2, the sign of the second-order gradients generated by the Laplacian operator can tell us whether the central pixel lies on the dark or light side of an edge if it is on the edge. Let G2N(x0, y0) and G2D(x0, y0) represent the second-order gradients generated by the Laplacian operator LN and LD (9), respectively. Let G2(x0, y0) be the summation of the absolute values of these two measures (23).

Fig. 9. Test image ‘‘Lena’’ contaminated by AWGN with radd = 10 and restored by the bilateral filter with different parameter selections.

820

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

G2N ðx0 ; y0 Þ ¼ LN  Xðx0 ; y0 Þ; G2D ðx0 ; y0 Þ ¼ LD  Xðx0 ; y0 Þ G2 ðx0 ; y0 Þ ¼ jG2N ðx0 ; y0 Þj þ jG2D ðx0 ; y0 Þj

ð23Þ

Our first criterion for a pixel to be claimed as an edge pixel is that two second-order gradients must agree with each other on the dark-or-light indication. In other words, their product must be positive (24).

G2N ðx; yÞ  G2D ðx; yÞ > 0

ð24Þ

As reported in Section 2, the image noise estimation process is performed only on the homogeneous areas of the noisy image. The homogeneous blocks are identified by the P% principle where the threshold value was based on the histogram of the homogeneity measure. Here, we borrow the same idea to identify the most notable edges for enhancement. It is assumed that Pen% of the whole image will be processed as edges and enhanced by the proposed algorithm. The histogram of G2 (23) is computed first. Then, the threshold value Gth 2 is selected to be the G2 value when the accumulated histogram reaches (100  Pen)% of the whole image. Let f(G2) be the discrete distribution of G2 and F G2 ðtÞ be the discrete cumulative distribution of G2. The threshold Gth 2 is defined by (25).

F G2 ðtÞ ¼

X

f ðG2 Þ

G2 6t

F G2 ðGth 2 Þ ¼



WH W 2b

 en Þ ; bcis the floor function:  ð100P 100

ð25Þ Fig. 10. The edge map of test image Lena (rn = 10) generated by the Laplacian operator and the rule defined in (28) with Pen = 10.

Our second criterion for a pixel to be claimed as an edge pixel is that its second-order gradient G2 is greater than or equal to the threshold value of Gth 2 . Since the noise estimation process is block-based, our edge searching process is also block-based. Let Bedge represents the subxy set of Bxy (2) which satisfies the first and second criteria. Our third criterion to pick up an edge pixel is that the cardinality of Bedge xy must be greater than one (26).

Bedge xy

8

9

G2N ðm; nÞ  G2D ðm; nÞ > 0 > >

< =

¼ ðm; nÞ G2 ðm; nÞ P Gth 2

> > : ;

ðm; nÞ 2 Bxy

ð26Þ

jBedge xy j > 1 Let Ben represents the subset of image pixels which are claimed to be edge pixels and will be enhanced by the noise reduction process; it can be defined by the following equation.

Ben

8

9

G2N ðx; yÞ  G2D ðx; yÞ > 0 > > <

=

G ðx; yÞ P Gth ¼ ðx; yÞ 2 2

edge > > : ;

jB j > 1

ð27Þ

xy

An edge map defined by (28) is shown in Fig. 10. As it can be seen, the strong edges are clearly identified by (27) and the weak edges are ruled out. Although there are still some isolated false edges, the enhancement on these pixels will be bound by the offset discussed in the following. Therefore, the enhancement process will not result in new artificial noise.

8 ðx; yÞ 2 Ben andG2N ðx; yÞ < 0 and G2D ðx; yÞ < 0 > < 16 Iðx; yÞ ¼ 127 ðx; yÞ R Ben > : 235 ðx; yÞ 2 Ben and G2N ðx; yÞ > 0 and G2D ðx; yÞ > 0 ð28Þ Let IMAX, IMIN, and IAVG represent the highest, lowest, and average intensity of the filtering window X(x0, y0)(16), respectively (29-31).

IMAX ¼ fMAXðIðx; yÞÞjðx; yÞ 2 Xðx0 ; y0 Þg

ð29Þ

IMIN ¼ fMINðIðx; yÞÞjðx; yÞ 2 Xðx0 ; y0 Þg

ð30Þ

IAVG ¼

1 Wf  Wf

X

Iðx; yÞ

ð31Þ

ðx;yÞ2Xðx0 ;y0 Þ

We propose to select the offset s(x0, y0) by the following rule where ren is a user-provided parameter. 8 ðx0 ; y0 Þ 2 Ben and G2N ðx0 ; y0 Þ > 0 > < MINðIMAX  IAVG ; ren rest Þ sðx0 ; y0 Þ ¼ MAXðIMIN  IAVG ; ren rest Þ ðx0 ; y0 Þ 2 Ben and G2N ðx0 ; y0 Þ < 0 > : 0 otherwise

ð32Þ The rationales of our algorithm (32) are explained as follows. First, the sharpening should be performed only when the central pixel is really on an edge, otherwise new noise will be introduced. Secondly, the sharpening should be performed with proper extent. If the central pixel was already on a sharp edge, too much enhancement will generate annoying result. Therefore, the sharpening should be bound by the estimated noise level. On the other hand, if the image was really noisy, the sharpening should be bound by the extent of the intensity variations inside the filter window. To show the performance of the proposed modified Bilateral filter, the test image Lena contaminated by AWGN with radd = 10 and restored by the bilateral filter and the modified bilateral filter with sharpening offset is shown in Fig. 11. As it can be seen from the difference image (Fig. 11(d)), the details of the hat as well as the edges around her eyes are enhanced. Another test image Barbara contaminated by AWGN with radd = 10 and restored by the bilateral filter and the modified bilateral filter with sharpening offset is shown in Fig. 12. Instead of displaying the whole image, two zoomed-in areas are presented to demonstrate the enhancement more clearly. Notice that these two filtered images achieve the same PSNR values, but the visual effects are quite different. From the two examples given here we can conclude that our proposed modified bilateral filter provide a very flexible solution for image noise reduction and edge enhancement. Users can obtain different filtered images by manipulating the combinations of four parameters, i.e. the filtering window size Wf, the ‘‘strength’’ of the range filter r, the enhancement percentage Pen, and the enhancement strength ren.

821

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

Fig. 11. Test image ‘‘Lena’’ contaminated by AWGN with radd = 10 and restored by the bilateral filter and the modified bilateral filter with sharpening offset.

6. Non-local mean filter Most of the well-known averaging filters can be classified as local averaging (or mean) filters. They are ‘‘local’’ because the averaging process is performed on pixels located in the neighborhood, usually a small window. Also, the weights of the averaging depend on the closeness or similarity between two pixels. Baudes, Coll and Morel [30–32] proposed a new averaging filter which is now known as the non-local mean filter (NLM). The non-local mean filter differs from the traditional local averaging filters in two ways. First, it performs averaging on pixels belonging to the so-called searching window XSH (Fig. 13). The window size of XSH is usually larger than the size of local averaging filters. What’s more important is that the weights of the averaging depend on the similarity between two blocks known as similarity windows (XSI) (Fig. 13). While the local averaging filters look for the similarity between gray levels of two pixels, the NLM filter looks for the similarity between structures of two blocks (or windows). 6.1. The filter formulation To show the resemblance between the NLM filter and the bilateral filter, we try to formulate the NLM filter by the same notations used in Eqs. (15)–(19). The NLM filter is defined mainly by (33).

g^NLM ðx0 ; y0 Þ ¼

X 1 C NLM ðx0 ; y0 Þ ðx;yÞ2X ðx SH

kNLM ðx; y; x0 ; y0 ÞIðx; yÞ

ð33Þ

0 ;y0 Þ

XSH(x0, y0) defined by (34) is the searching window centered at (x0,y0).

XSH ðx0 ; y0 Þ ¼ fðx; yÞjjx  x0 j 6 W hSH and jy  y0 j 6 W hSH g

ð34Þ

WSH is the window size of the searching window and WhSH is half of the size.

W SH ¼ 2W hSH þ 1

ð35Þ

The weight kNLM(x, y; x0, y0) is a Gaussian kernel which exploits the similarity between two similarity windows XSI(x0, y0) and XSI(x, y).

 v ðx; y; x0 ; y0 Þ kNLM ðx; y; x0 ; y0 Þ ¼ exp  2 h X v ðx; y; x0 ; y0 Þ ¼ ðx1 ; y1 Þ 2 XSI ðx0 ; y0 Þ ðx2 ; y2 Þ 2 XSI ðx; yÞ

ð36Þ

jx1  x2 j2 þ jy1  y2 j2 jIðx1 ; y1 Þ  Iðx2 ; y2 Þj2  exp  2q2

!

ð37Þ The similarity window is defined by (38) and (39) where WSI is the window size of the similarity window and WhSI is half of the size.

822

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

Fig. 12. Test image ‘‘Barbara’’ contaminated by AWGN with radd = 10 and restored by the bilateral filter and the modified bilateral filter with sharpening offset.

XSI ðx0 ; y0 Þ ¼ fðx; yÞjjx  x0 j 6 W hSI and jy  y0 j 6 W hSI g

ð38Þ

W SI ¼ 2W hSI þ 1

ð39Þ

CNLM(x0, y0) defined by (40) is a normalization factor.

C NLM ðx0 ; y0 Þ ¼

X

kNLM ðx; y; x0 ; y0 Þ

ð40Þ

ðx;yÞ2XSH ðx0 ;y0 Þ

6.2. Searching for the best parameter settings There are four parameters (i.e. WSH, WSI, h and q) which decide the performance of the NLM filter. We did the same exhaustive searching experiments as we did for the bilateral filter. The following rules of thumb are observed. 1. If WSH and WSI are fixed, while (h, q) = (4, 1.5) achieves the highest PSNR at low noise levels, (4, 2) or (6, 2) achieves the highest PSNR at medium to high noise levels (Table 1). However, the differences are negligible.

2. When the other three parameters were fixed, larger WSI resulted in higher PSNR. But again, the improvements are negligible in terms of PSNR and visual quality (Table 2). Therefore, WSI can be fixed at 7 and q can be set to 1.5. 3. h can easily result in blurring if it is set to a large value. 4 is a good compromise. 4. The most important parameter which has great impact to both the visual quality and computation complexity is WSH. It can be set to a value between 7 and 11 if the noise level is low, and a value between 13 and 15 if the noise level is medium (Table 3). For highly noisy images, it will take 21 to achieve acceptable visual quality. In a short conclusion, WSH can be decided by (41).

W SH

8 7 > > < ¼ 9 > > : the odd integer greater than and closest to

rn 6 7 7 < rn 6 10

rest 10 < rn ð41Þ

823

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

6.3. Reducing the Complexity The complexity of the NLM filter is order of (NImg, WSH, WSI) where NImg is the image size. Although the NLM filter has state of the art noise filtering performance, the extremely high complexity makes it impossible for real application. Many algorithms have been proposed to reduce its complexity, for example [33,34].

One thing we learned from the experiments is that the searching window plays an important role on the visual quality of noise reduction. We also believe that searching in the right direction is more important than searching in a larger window. As mentioned in Section 3.1, the angle of the gradient vector produced by the Sobel operator tells us the direction of an edge. Let’s produce an edge map by the following definitions where axy is defined by (6)-(8).

Table 2 Performance comparisons of the original NLM(BCM) with three parameters fixed (q and h and WSH) and only one parameter varied (WSI). Noisy image Lena Lena Lena Lena Lena Lena

with with with with with with

rn = 10 rn = 10 rn = 10 rn = 20 rn = 20 rn = 20

WSI

WSH

q

h

PSNR

5 7 9 5 7 9

9 9 9 19 19 19

1.5 1.5 1.5 2 2 2

4 4 4 4 4 4

35.95 36.09 36.12 32.02 32.31 32.44

Table 3 Performance comparisons of the original NLM (BCM) with three parameters fixed (q and h and WSI) and only one parameter varied (WSH). Noisy image

Fig. 13. Scheme of the non-local mean filter.

Table 1 Performance comparisons of the original NLM(BCM) with two parameters fixed (WSI and WSH) and two parameters varied (q and h). WSI

WSH

q

h

PSNR

Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena

with with with with with with with with with with with with with with with

rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5

5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

7 7 7 7 7 7 7 7 7 7 7 7 7 7 7

1 1 1 1 1 1.5 1.5 1.5 1.5 1.5 2 2 2 2 2

2 4 6 8 10 2 4 6 8 10 2 4 6 8 10

38.34 39.05 38.11 37.07 36.14 37.75 38.98 38.5 37.72 36.94 37.36 38.79 38.55 37.93 37.24

Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena

with with with with with with with with with with with with with with with

rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15 rn = 15

5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

7 7 7 7 7 7 7 7 7 7 7 7 7 7 7

1 1 1 1 1 1.5 1.5 1.5 1.5 1.5 2 2 2 2 2

2 4 6 8 10 2 4 6 8 10 2 4 6 8 10

30.84 33.49 32.75 31.81 31.06 30.21 33.74 33.36 32.5 31.74 29.76 33.68 33.6 32.82 32.07

Noisy image

Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena Lena

with with with with with with with with with with with with with with with with

rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 5 rn = 20 rn = 20 rn = 20 rn = 20 rn = 20 rn = 20 rn = 20 rn = 20

WSI

WSH

q

h

PSNR

7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7

7 9 11 13 15 17 19 21 7 9 11 13 15 17 19 21

1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 2 2 2 2 2 2 2 2

4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

39.01 39.01 38.98 38.95 38.9 38.87 38.83 38.8 32.08 32.39 32.5 32.51 32.46 32.39 32.31 32.23

Fig. 14. The edge map generated by the Sobel operator and the rule defined in (42).

824

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

(x0,y0) (xx0,yy0)

( 0,yy0) (x

central quarter

1st 1 quadrant

4th quadrant

Fig. 15. Searching windows for our proposed QBCM filter.

Table 4 Performance comparisons between the original NLM (BCM) and our proposed algorithm (QBCM) in terms of PSNR. Noisy image

WSI

WSH

q

h

QBCM

BCM

Mandrill with rn = 5 Mandrill with rn = 10 Mandrill with rn = 15 Mandrill with rn = 20 Lena with rn = 5 Lena with rn = 10 Lena with rn = 15 Lena with rn = 20

7 7 7 7 7 7 7 7

7 9 15 21 7 9 15 21

1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5

4 4 4 4 4 4 4 4

31.96 30.09 28.10 26.27 38.12 35.18 33.04 31.33

33.04 30.66 28.06 26.10 38.99 36.08 33.72 31.85

8 16 > > > < 235 Iðx; yÞ ¼ > 90 > > : 160

jaxy j 6 100 jaxy j P 800 100 < axy < 800 800 < axy < 100

ð42Þ

The ‘‘edge’’ map of the test image ‘‘Cameraman’’ corrupted by AWGN with rn = 10 is shown in Fig. 14. Let’s focus on the tripod. It can be seen that the angle of the Sobel gradient really can catch the directions of the edges shown with different gray levels. Therefore, we can choose the Sobel operator as our homogeneity analyzer and define the searching area based on the following rule. The rationale of this rule is very straightforward.

Fig. 16. Restored images of ‘‘Mandrill’’ contaminated by AWGN with radd = 5 and the difference image.

825

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

Fig. 17. Restored images of ‘‘Lena’’ contaminated by AWGN with radd = 20 and the difference image.

Table 5 Performances of the modified bilateral filter and the original bilateral filter applied to the same test image with different resolutions. Noisy image

Wf

rd

r

Pen

ren

Modified bilateral

Bilateral

Lena (512  512) with rn = 10 Lena (256  256) with rn = 10

9 9

4 4

3 3

20 20

2 2

32.71 31.20

32.97 30.78

Table 6 Performances of QBCM and BCM applied to the same test image with different resolutions. Noisy image

WSI

WSH

q

h

QBCM

BCM

Mandrill (512  512) with rn = 5 Mandrill (256  256) with rn = 5

7 7

7 9

1.5 1.5

4 4

31.96 26.97

33.04 27.59

8 100 < axy < 800 > < 1st Quadrant XqSH ðx0 ; y0 Þ ¼ 4th Quadrant 800 < axy < 100 > : central quarter jaxy j 6 100 orjaxy j P 800

ð43Þ

The terms of ‘‘1st Quadrant’’, ‘‘4th Quadrant’’, and ‘‘central quarter’’ are referred to one fourth of the original searching window and illustrated in Fig. 15. Based on this quartered rule, the complexity of the original NLM filter is reduced by 75%. We will refer to the original NLM filter by Baudes, Coll and Morel as the BCM, and our proposed modified NLM filter with only 25% complexity as the QBCM.

The performance of our proposed algorithm in terms of PSNR can be shown by the data given in Table 4. The performance in terms of visual quality can be shown by examples given in Figs. 16 and 17. For images with fine details like ‘‘Mandrill’’, the QBCM loses about 1 dB at low noise level (rn = 5) (Table 4). Despite of the defeat, QBCM achieves the same visual quality as BCM does (Fig. 16 (a) and (b)). The difference only lies on the regions where there are even no weak edges (Fig. 16(c)) and Sobel edge detectors can’t help. Fortunately, human eyes are not sensitive to details in these regions. For medium to high noise levels, QBCM achieves almost the same PSNRs as the BCM does and sometimes even better. This particular example tells us that more samples do not always provide better denoising results. Sometimes, bad samples may degrade the restoration. For image ‘‘Lena’’, the difference between BCM and QBCM remains in 1dB. The difference image shown in Fig. 17(c) tells us that BCM and QBCM provide the same noise reduction performance in terms of visual quality. All these observation confirms our theory

826

S.-M. Yang, S.-C. Tai / J. Vis. Commun. Image R. 23 (2012) 812–826

that searching in the right directions is far more important than searching in a larger area. 6.4. Comments on image resolution Both the bilateral filter and the non-local mean filter are spatialdomain filters. They are also averaging filters which replace the central pixel with the average of its neighbors to reduce image noise. The success of noise reduction by averaging relies on the correlation between the central pixel and its neighbors. When a picture is taken with higher resolution, the correlation between pixels also increases. Therefore, we expect to see better performance when an averaging filter is applied to the same test image with higher resolution. As shown in Tables 5 and 6, both the bilateral and the non-local mean filter achieve higher PSNR when applied to images with higher resolution. For images with fine details such as Mandrill, the gain on PSNR due to higher resolution is even more notable. This is important because the resolution of modern digital cameras is getting higher and higher. 7. Conclusions We proposed a design framework for hybrid approaches of image noise estimation and reduction, assuming that it is additive white Gaussian noise (AWGN). Four high-pass operators were evaluated for being the homogeneity analyzer in the estimation process. Their performances on noise estimation were compared by Monte–Carlo simulations conducted on a variety of test images over a large range of noise levels. The experimental results showed that all of them make good noise estimators. Two well-known nonlinear filters, the bilateral and the non-local mean, were reviewed. How to integrate the noise estimation with noise reduction was then investigated. The best parameter settings of these two filters were searched thoroughly by simulations. The impact of each parameter to the filtering performance was discussed. A modified bilateral filter with edge enhancement was proposed. The enhancement was achieved by incorporating side information from the noise estimator. A modified non-local mean filter with lower complexity was also present. By searching in the right direction which was pointed out by the noise estimator, the complexity of noise filtering was reduced tremendously but the performance was maintained or even improved. References [1] J.-S. Lee, Digital image smoothing and the sigma filter, Computer Vision, Graphics, and Image Processing 24 (2) (1983) 255–269. [2] G.A. Mastin, Adaptive filters for digital image noise smoothing: An evaluation, Computer Vision, Graphics, and Image Processing 31 (1) (1985) 103–121. [3] C. Tomasi, R. Manduchi, Bilateral filtering for gray and color images, in: Sixth International Conference on Computer Vision, New Delhi, India, 1998, pp. 839– 846. [4] F. Russo, A method based on piecewise linear models for accurate restoration of images corrupted by Gaussian noise, IEEE Transactions on Instrumentation & Measurement 55 (6) (2006) 1935–1943. [5] F. Russo, Technique for image denoising based on adaptive piecewise linear filters and automatic parameter tuning, IEEE Transactions on Instrumentation & Measurement 55 (6) (2006) 1362–1367.

[6] F. Russo, An image-enhancement system based on noise estimation, in IEEE Transactions on Instrumentation & Measurement 56 (4) (2007) 1435–1442. [7] B. Zhang, J.P. Allebach, Adaptive bilateral filter for sharpness enhancement and noise removal, IEEE Transactions on Image Processing 17 (5) (2008) 664–678. [8] M. Ghazal, A. Amer, A. Ghrayeb, Structure-oriented multidirectional wiener filter for denoising of image and video signals, Circuits and Systems for Video Technology, IEEE Transactions on 18 (12) (2008) 1797–1802. [9] W.T. Freeman, E.C. Pasztor, O.T. Carmichael, Learning low-level vision, International Journal of Computer Vision 40 (1) (2000) 25–47. [10] D. Lowe, Object recognition from local scale invariant features, in: Proceedings of the IEEE International Conference Computer Vision, 1999, pp. 1150–1157. [11] G. Calvagno, F. Fantozzi, R. Rinaldo, A. Viareggio, Model-based global and local motion estimation for videoconference sequences, IEEE Transactions on Circuits and Systems for Video Technology 14 (9) (2004) 1156–1161. [12] S.I. Olsen, Estimation of noise in images: An evaluation, CVGIP: Graphical Models and Image Processing 55 (4) (1993) 319–323. [13] D.L. Donoho, J.M. Johnstone, Ideal spatial adaptation by wavelet shrinkage, Biometrika 81 (3) (1994) 425–455. [14] H.H. Khalil, R.O.K. Rahmat, W.A. Mahmoud, Estimation of noise in gray-scale and colored images using median absolute deviation (MAD), in: Proceedings of the 2008 3rd International Conference on Geometric Modeling and Imaging, 2008, pp. 92–97. [15] M. Salmeri, A. Mencattini, E. Ricci, A. Salsano, Noise estimation in digital images using fuzzy processing, in: Proceedings of 2001 International Conference on Image Processing, 2001, pp. 517–520. [16] J. Immerkær, Fast noise variance estimation, Computer Vision and Image Understanding 64 (2) (1996) 300–302. [17] K. Rank, M. Lendl, R. Unbehauen, Estimation of image noise variance, IEE Proceedings on Vision, Image and Signal Processing 146 (2) (1999) 80–84. [18] R.C. Bilcu, M. Vehvilainen, New method for noise estimation in images, in: Proceedings of IEEE-Eurasip International Workshop on Nonlinear Signal and Image Processing, 2005, pp. 290–293. [19] A. Amer, E. Dubois, Fast and reliable structure-oriented video noise estimation, IEEE Transactions on Circuits & Systems for Video Technology 15 (1) (2005) 113–118. [20] S. Aja-Fernandez, G. Vegas-Sanchez-Ferrero, M. Martin-Fernandez, C. AlberolaLopez, Automatic noise estimation in images using local statistics. Additive and multiplicative cases, Image and Vision Computing 27 (6) (2009) 756– 770. [21] D.-H. Shin, R.-H. Park, S. Yang, J.-H. Jung, Block-based noise estimation using adaptive gaussian filtering, IEEE Transactions on Consumer Electronics 51 (1) (2005) 218–226. [22] S.M. Yang, S.C. Tai, Fast and reliable image noise estimation using a hybrid approach, Journal of Electronic Imaging 19 (3) (2010). 033007-1. [23] S.M. Yang, S.C. Tai, A Fast and Reliable Algorithm for Video Noise Estimation Based on Spatio-Temporal Sobel Gradients, in: the 1st International Conference on Electrical, Control and Computer Engineering, ECCE2011, June 21–22, 2011, Pahang, Malaysia. [24] R.C. Gonzalez, R.E. Woods, Digital Image Processing, 2nd ed., Prentice-Hall, Inc., New Jersey, USA, 2002. [25] S.R. Gunn, On the discrete representation of the Laplacian of Gaussian, Pattern Recognition 32 (1999) 1463–1472. [26] S. Paris, P. Kornprobst, J. Tumblin, F. Durand, A gentle introduction to bilateral filtering and its application, presented at the ACM SIGGRAPH, 2007. [27] M. Zhang, B. Gunturk. A new image denoising method based on the bilateral filter, in: IEEE International Conference on Acoustics, Speech, Signal Processing (ICASSP), 2008, pp. 9299–9232. [28] M. Zhang, B.K. Gunturk, Multiresolution bilateral filtering for image denoising, IEEE Transactions on Image Processing 17 (12) (2008) 2324–2333. [29] B.K. Gunturk, Bilateral filter: theory and application, in: R. Lukac (Ed.), Computational Photography, CRC Press, 2011. [30] A. Baudes, B. Coll, J.M. Morel, A review of image denoising algorithms, with a new one, Multiscale Modeling & Simulation (4) (2005) 490–530. [31] A. Buades, B. Coll, J.M. Morel, A non-local algorithm for image denoising, in: IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR), 2005, pp. 60–65. [32] A. Buades, B. Coll, J.M. Morel, Nonlocal image and movie denoising, International Journal of Computer Vision 76 (2) (2008) 123–139. [33] M. Mahmoudi, G. Sapiro, Fast image and video denoising via nonlocal means of similar neighborhoods, IEEE Signal Processing Letters 12 (2) (2005) 839–842. [34] R.C. Bilcu, M. Vehvilainen, Modified non-local averaging for image de-noising, WSEAS Transactions on Circuits and Systems 5 (7) (2006).