计算机科学与技术（参考文献）A Robust and Fast Non-local Algorithm for Image Denoising

Liu YL,Wang J,Chen X et al.A robust and fast non-local means algorithm for image denoising.JOURNAL OF COM- PUTER SCIENCE AND TECHNOLOGY 23(2):270-279 Mar.2008 A Robust and Fast Non-Local Means Algorithm for Image Denoising Yan-Li Liu.2(刘艳丽)，Jin Wang(王进)，Xi Chen(陈曦)，Yan-Wen Guo3(郭延文) and Qun-Sheng Peng,2(彭群生) State Key Lab of CAD&CG,Zhejiang University,Hangzhou 310058,China 2Department of Mathematics,Zhejiang University,Hangzhou 310058,China 3State Key Lab of Novel Software Technology,Nanjing University,Nanjing 210000,China E-mail:liuyanli@cad.zju.edu.cn Revised January 7,2008 Abstract In the paper,we propose a robust and fast image denoising method.The approach integrates both Non- Local means algorithm and Laplacian Pyramid.Given an image to be denoised,we first decompose it into Laplacian pyramid.Exploiting the redundancy property of Laplacian pyramid,we then perform non-local means on every level image of Laplacian pyramid.Essentially,we use the similarity of image features in Laplacian pyramid to act as weight to denoise image.Since the features extracted in Laplacian pyramid are localized in spatial position and scale,they are much more able to describe image,and computing the similarity between them is more reasonable and more robust.Also,based on the efficient Summed Square Image (SSI)scheme and Fast Fourier Transform (FFT),we present an accelerating algorithm to break the bottleneck of non-local means algorithm-similarity computation of compare windows.After speedup,our algorithm is fifty times faster than original non-local means algorithm.Experiments demonstrated the effectiveness of our algorithm. Keywords image denoising,non-local means,Laplacian pyramid,summed square image,FFT 1 Introduction ods.Whereas,generally speaking,the information en- coded in a natural image is redundant to some extent, With the prevalence of digital cameras and scan- that is to say,there may exist some repeat patterns ners,digital images can be easily acquired nowadays. in natural images,particularly in textured or periodic Unfortunately,digital images are often degraded by case.Based on this observation,Buadeslio]developed noise during acquisition or transmission process.Var- a non-local image denoising algorithm which takes full ious image-related applications,such as medical im- advantage of image redundancy.For the sake of sim- age analysis,image segmentation,and object detec- plicity,we name it the spatial nl-means in the paper. tion,etc.,generally require effective noise suppression Like many noise reduction algorithm,the method is method to further produce reliable results.Therefore, also based on weighted average.The essence of the image denoising has been one of the most important method is:to estimate a certain pixel,the method uses and widely studied problems in image processing and the similarities between it and all the other pixels in computer vision.So far,image denoising methods image to act as weight,and the similarities are not can be basically divided into two categories:spatial computed from pixels themselves but from their neigh- filtering methods(1-4 and transform domain filtering boring window (compare window).The algorithm has methods5-9] demonstrated strong superiority over local-based spa- In spatial domain,for every pixel of noisy image, tial methods such as Gaussian filter,bilateral filter in many existing noise reduction methods employ the in- terms of both PSNR and visual quality.However,cer- formation of its nearby local region to estimate its de- tain critical issues need further investigation.Firstly, noised version.Examples include Gaussian filterl],me- the computational cost is high for its pixel-by-pixel win- dian filterl3),bilateral filterl4)and so on.We call this dow matching,and it limits the method to be widely kind of denoising algorithms local-based spatial meth- used in application.Secondly,the quality and compu- Regular Paper This work is supported by the National Grand Fundamental Research 973 Program of China (Grant No.2002CB312101),the National Natural Science Foundation of China (Grant Nos.60403038 and 60703084)and the Natural Science Foundation of Jiangsu Province (Grant No.BK2007571)

Liu YL, Wang J, Chen X et al. A robust and fast non-local means algorithm for image denoising. JOURNAL OF COM￾PUTER SCIENCE AND TECHNOLOGY 23(2): 270–279 Mar. 2008 A Robust and Fast Non-Local Means Algorithm for Image Denoising Yan-Li Liu1,2 (刘艳丽), Jin Wang1 (王 进), Xi Chen1 (陈 曦), Yan-Wen Guo3 (郭延文) and Qun-Sheng Peng1,2 (彭群生) 1State Key Lab of CAD & CG, Zhejiang University, Hangzhou 310058, China 2Department of Mathematics, Zhejiang University, Hangzhou 310058, China 3State Key Lab of Novel Software Technology, Nanjing University, Nanjing 210000, China E-mail: liuyanli@cad.zju.edu.cn Revised January 7, 2008. Abstract In the paper, we propose a robust and fast image denoising method. The approach integrates both Non￾Local means algorithm and Laplacian Pyramid. Given an image to be denoised, we first decompose it into Laplacian pyramid. Exploiting the redundancy property of Laplacian pyramid, we then perform non-local means on every level image of Laplacian pyramid. Essentially, we use the similarity of image features in Laplacian pyramid to act as weight to denoise image. Since the features extracted in Laplacian pyramid are localized in spatial position and scale, they are much more able to describe image, and computing the similarity between them is more reasonable and more robust. Also, based on the efficient Summed Square Image (SSI) scheme and Fast Fourier Transform (FFT), we present an accelerating algorithm to break the bottleneck of non-local means algorithm — similarity computation of compare windows. After speedup, our algorithm is fifty times faster than original non-local means algorithm. Experiments demonstrated the effectiveness of our algorithm. Keywords image denoising, non-local means, Laplacian pyramid, summed square image, FFT 1 Introduction With the prevalence of digital cameras and scan￾ners, digital images can be easily acquired nowadays. Unfortunately, digital images are often degraded by noise during acquisition or transmission process. Var￾ious image-related applications, such as medical im￾age analysis, image segmentation, and object detec￾tion, etc., generally require effective noise suppression method to further produce reliable results. Therefore, image denoising has been one of the most important and widely studied problems in image processing and computer vision. So far, image denoising methods can be basically divided into two categories: spatial filtering methods[1−4] and transform domain filtering methods[5−9] . In spatial domain, for every pixel of noisy image, many existing noise reduction methods employ the in￾formation of its nearby local region to estimate its de￾noised version. Examples include Gaussian filter[1], me￾dian filter[3], bilateral filter[4] and so on. We call this kind of denoising algorithms local-based spatial meth￾ods. Whereas, generally speaking, the information en￾coded in a natural image is redundant to some extent, that is to say, there may exist some repeat patterns in natural images, particularly in textured or periodic case. Based on this observation, Buades[10] developed a non-local image denoising algorithm which takes full advantage of image redundancy. For the sake of sim￾plicity, we name it the spatial nl-means in the paper. Like many noise reduction algorithm, the method is also based on weighted average. The essence of the method is: to estimate a certain pixel, the method uses the similarities between it and all the other pixels in image to act as weight, and the similarities are not computed from pixels themselves but from their neigh￾boring window (compare window). The algorithm has demonstrated strong superiority over local-based spa￾tial methods such as Gaussian filter, bilateral filter in terms of both PSNR and visual quality. However, cer￾tain critical issues need further investigation. Firstly, the computational cost is high for its pixel-by-pixel win￾dow matching, and it limits the method to be widely used in application. Secondly, the quality and compu￾Regular Paper This work is supported by the National Grand Fundamental Research 973 Program of China (Grant No. 2002CB312101), the National Natural Science Foundation of China (Grant Nos. 60403038 and 60703084) and the Natural Science Foundation of Jiangsu Province (Grant No. BK2007571)

Yan-Li Liu et al:A Robust and Fast Non-Local Means Algorithm for Image Denoising 271 tation cost of the algorithm is closely related to the size Image(SSI)scheme and Fast Fourier Transform(FFT). of compare window which represents geometrical con- the per-pixel neighborhood matching in Laplacian im- figuration of pixel's neighboring region.For small com- age is converted into the SSI pre-computing and ef- pare window,the algorithm is restricted to suppresing ficient FFT.Computational complexity analysis and the high-frequency noise and cannot remove the low- experiments indicate that our accelerated algorithm is frequency (large-scale)noise.For large compare win- more than 50 times faster than the original non-local dow,the algorithm removes the low-frequency noise ef algorithm.Using our algorithm,it normally takes less fectively,but becomes less sensitive to small details and than 1 second to denoise a 640 x 480 image,and less the additional burden on the computational resources than 11 seconds to denoise a 500 Megabyte photograph, becomes unacceptable. which enables the algorithm applicable to practical sit- Since the essential goal of denoising is to preserve uations. the image features while reducing noise effectively,a The main contributions of our paper are as follows. logical extension of spatial denoising is to transform 1)We propose a robust non-local denoising method images into a representation that distinctive features based on Laplacian pyramid.Combining Laplacian such as edges can be extracted from image and per- pyramid and non-local means for image denoising has form noise reduction algorithms in this domain.As for the following advantages.(i)Laplacian pyramid decom- denoising in frequency domain,we are more familiar pose noisy image into easily manipulatable image com- with waveletsl5-9]which fall into two classes:orthog- ponents.By performing non-local means (nl-means)on onal wavelet and non-orthogonal wavelet.Orthogonal different levels of Laplacian pyramid with different sizes wavelet compositions are critically sampled (i.e.,nonre- of compare windows,we not only reduce high-frequency dundant)image descriptions.It has been successfully noise but also reduce low-frequency noise while preserv- used in image compression1,121 and threshold-based ing the image details (edges,textures)well.(ii)Since image denoising5,131 However,the non-redundant the features extracted in Laplacian pyramid are local- property of orthogonal wavelet means they are not suit- ized in spatial position and scale,they are much more able for non-local means denoising which builds upon able to describe image,and computing the similarity image redundancy. between them is more reasonable and more robust.(iii) In this paper,we emphasis on overcomplete non- As the level of Laplacian pyramid increases,the amount orthogonal representations of image. Some non- of noise in Laplacian image decreases rapidly,therefore, orthogonal wavelet representations are redundant to denoising in Laplacian image is less noise-sensitive than some extent.It should be noted that wavelet decom- in spatial domain. poses image into three sub-bands for every level,im- 2)We also present an accelerating algorithm to plying that connected spatial structures analyzed at in speed up similarity computation of compare windows. creasing resolutions are split into separate sub-bands. Based on efficient SSI scheme and FFT.the accelerated thereby losing their spatial connectivity4.Another algorithm is fifty times faster than the spatial nl-means drawback of wavelet representation is that comparing The result of spatial nl-means leaves a trail of smear. the similarity respectively in three sub-bands will in- When coming across heavy noise,the spatial nl-means duce heavy computational cost.Based on these con- introduces artifacts such as canvas effect which is pro- siderations,Laplacian pyramid5],being a redundant nounced in flat regions and soft edges of denoised image. image representation,is selected in the paper.Due to In contrast,our method reduces noise effectively and the overcompletion property of Laplacian image,every preserves the image details and weak edges well.As frequency component image is redundant,thus we can noise level increases,the difference of results between compute the non-local similarity in frequency compo- our method and the spatial nl-means is even more strik- nent.In the paper,we first exploit Laplacian pyra- ing. mid to decompose the given noisy image and generate The rest of this paper is organized as follows.Section a series of manipulatable components,which allows us 2 gives background and related work about nl-means to develop different stategies on different components. algorithm.Laplacian pyramid and how our algorithm Then we perform non-local means on every level of integrates both nl-means algorithm and Laplacian pyra- Laplacian pyramid with different sizes of compare win- mid as well as our acceleration scheme are detailed in dows. Section 3.Experimental results are demonstrated in Also,an acceleration method to our algorithm is pro- Section 4.Section 5 gives the conclusion of the whole posed in the paper.Exploiting the Summed Squares paper

Yan-Li Liu et al.: A Robust and Fast Non-Local Means Algorithm for Image Denoising 271 tation cost of the algorithm is closely related to the size of compare window which represents geometrical con- figuration of pixel’s neighboring region. For small com￾pare window, the algorithm is restricted to suppresing the high-frequency noise and cannot remove the low￾frequency (large-scale) noise. For large compare win￾dow, the algorithm removes the low-frequency noise ef￾fectively, but becomes less sensitive to small details and the additional burden on the computational resources becomes unacceptable. Since the essential goal of denoising is to preserve the image features while reducing noise effectively, a logical extension of spatial denoising is to transform images into a representation that distinctive features such as edges can be extracted from image and per￾form noise reduction algorithms in this domain. As for denoising in frequency domain, we are more familiar with wavelets[5−9] which fall into two classes: orthog￾onal wavelet and non-orthogonal wavelet. Orthogonal wavelet compositions are critically sampled (i.e., nonre￾dundant) image descriptions. It has been successfully used in image compression[11,12] and threshold-based image denoising[5,13]. However, the non-redundant property of orthogonal wavelet means they are not suit￾able for non-local means denoising which builds upon image redundancy. In this paper, we emphasis on overcomplete non￾orthogonal representations of image. Some non￾orthogonal wavelet representations are redundant to some extent. It should be noted that wavelet decom￾poses image into three sub-bands for every level, im￾plying that connected spatial structures analyzed at in￾creasing resolutions are split into separate sub-bands, thereby losing their spatial connectivity[14]. Another drawback of wavelet representation is that comparing the similarity respectively in three sub-bands will in￾duce heavy computational cost. Based on these con￾siderations, Laplacian pyramid[15], being a redundant image representation, is selected in the paper. Due to the overcompletion property of Laplacian image, every frequency component image is redundant, thus we can compute the non-local similarity in frequency compo￾nent. In the paper, we first exploit Laplacian pyra￾mid to decompose the given noisy image and generate a series of manipulatable components, which allows us to develop different stategies on different components. Then we perform non-local means on every level of Laplacian pyramid with different sizes of compare win￾dows. Also, an acceleration method to our algorithm is pro￾posed in the paper. Exploiting the Summed Squares Image (SSI) scheme and Fast Fourier Transform (FFT), the per-pixel neighborhood matching in Laplacian im￾age is converted into the SSI pre-computing and ef- ficient FFT. Computational complexity analysis and experiments indicate that our accelerated algorithm is more than 50 times faster than the original non-local algorithm. Using our algorithm, it normally takes less than 1 second to denoise a 640 × 480 image, and less than 11 seconds to denoise a 500 Megabyte photograph, which enables the algorithm applicable to practical sit￾uations. The main contributions of our paper are as follows. 1) We propose a robust non-local denoising method based on Laplacian pyramid. Combining Laplacian pyramid and non-local means for image denoising has the following advantages. (i) Laplacian pyramid decom￾pose noisy image into easily manipulatable image com￾ponents. By performing non-local means (nl-means) on different levels of Laplacian pyramid with different sizes of compare windows, we not only reduce high-frequency noise but also reduce low-frequency noise while preserv￾ing the image details (edges, textures) well. (ii) Since the features extracted in Laplacian pyramid are local￾ized in spatial position and scale, they are much more able to describe image, and computing the similarity between them is more reasonable and more robust. (iii) As the level of Laplacian pyramid increases, the amount of noise in Laplacian image decreases rapidly, therefore, denoising in Laplacian image is less noise-sensitive than in spatial domain. 2) We also present an accelerating algorithm to speed up similarity computation of compare windows. Based on efficient SSI scheme and FFT, the accelerated algorithm is fifty times faster than the spatial nl-means. The result of spatial nl-means leaves a trail of smear. When coming across heavy noise, the spatial nl-means introduces artifacts such as canvas effect which is pro￾nounced in flat regions and soft edges of denoised image. In contrast, our method reduces noise effectively and preserves the image details and weak edges well. As noise level increases, the difference of results between our method and the spatial nl-means is even more strik￾ing. The rest of this paper is organized as follows. Section 2 gives background and related work about nl-means algorithm. Laplacian pyramid and how our algorithm integrates both nl-means algorithm and Laplacian pyra￾mid as well as our acceleration scheme are detailed in Section 3. Experimental results are demonstrated in Section 4. Section 5 gives the conclusion of the whole paper

272 J.Comput.Sci.Technol..Mar.2008,Vol.23,No.2 2 Background and Related Work In order to accelerate the algorithm.they introduced filters that eliminate unrelated neighborhoods from 2.1 Background weighted average.These filters are based on local aver- age gray values and gradients,preclassifying neighbor- In spatial nl-means,when modifying a pixel,the algorithm first computes the similarity between the hoods and reducing the influence of less-related areas when denoising a given pixel.Note that the times of ac- window centered around it and the windows centered celeration of Mahmoudi's algorithm is related to image around the other pixels in the whole image,then it takes the similarities as weight to adjust this pixel.It also can content,for example,it is 7 times faster than original spatial nl-means for some images,whereas,for some be interpreted that the method defines as"neighbor- other images it is more than 20 times faster than the hood of a pixel i"any set of pixels j in the image such original nl-means.In our method,after speedup,the that a window around j looks like a window around i, method is fifty times faster than the original nl-means and all pixels in that neighborhood are used for pre- for all the images. dicting the value at i.In essence,the method uses two windows'similarity other than two pixels'to estimate pixels and these windows slide on the whole image. 3 Our Robust and Fast Laplacian Pyramid If N2 is the number of pixels in the image,M2 is Based nl-Means the size of compare window,the complexity of this al- 3.1 Laplacian Pyramid gorithm is M2x N4.For computational purposes,the simplified algorithm restricts the search of similar win- Gaussian-Laplacian pyramid has been widely used dows in a“search window"of size S×S pixels.By in texture classifications.1,feature detection20l and fixing the search window at the size of 21 x 21 pixels so on.Given an image I(x,y),the Gaussian pyramid and the neighborhood size 7 x 7,the final complexity of of I is computed as follows: the algorithm is 49 x 441 x N2.However,even the sim plified algorithm still takes about 1 minute to denoise a Go(r,y）=I, 640 x 480 image on a common PC.It can be seen that Gk+1(x,y)=REDUCE(Gk(,y)),k=0,1,...,N. the high computational complexity makes it unfeasible (1) to tackle with practical issues. The REDUCE operation is carried out by convolving the image with a Gaussian low pass filter and then sub- sampling the filtered version at a factor of two.Repeat- 2.2 Related Work ing the REDUCE operation generates Gaussian pyra- To our knowledge,up to now,there is not much work mid.The filter mask is designed such that the center about nl-means denoising.Lately,Lukin proposed a pixel gets more weight than the neighboring ones and multiresolution-based approach to improve spatial nl- the remaining terms are chosen so that their sum is 1. means algorithm16l.The algorithm first performs nl- The separable Gaussian kernel is given by: means on noisy image in spatial domain with different w(r;c)=w(r)w(c). (2) block sizes and search ranges.Afterward,it transforms the resulting images (generally two images)into fre- As in the general case,in the paper,the kernel we used quency domain,and then combines the coefficients by is of size5×5,andw(r）=(1/16,1/4,3/8,1/4,1/16). another filter bank,finally inverses the filter bank to get Thus,the variance of the Gaussian kernel approximates the denoising result.It can be seen that this method 1.0. only considers getting rid of different frequency noise, In order to generate Laplacian pyramid,the second but it does not consider the redundancy of image fea- level of Gaussian image Gi is up-sampled by inserting tures itself and comparing similarity in transform do zeros after each pixel,smoothed with the small kernel as main.Conversely,we compute similarity between im- defined in Gaussian pyramid.These interpolation and age features at different scales,taking advantage of re filtering process is referred to as EXPAND operation. dundancy of different spatial frequency components of The difference between the previous level of Gaussian image.The methods can not only reduce different fre- image Go and expanded image Gi is named as Lapla- quency noise efficiently but also produce robust results. cian image.This process is continued to obtain a set of That is where our approach diverges from Lukin's al- band-pass filtered images.Expressed by formulation, gorithm. An algorithm for accelerating spatial nl-means de- Lk(z,y)=Gk(,y)-EXPAND(Gk+1(,y)), noising was proposed by Mahmoudi and Sapirol171. k=0,1,..,K (3)

272 J. Comput. Sci. & Technol., Mar. 2008, Vol.23, No.2 2 Background and Related Work 2.1 Background In spatial nl-means, when modifying a pixel, the algorithm first computes the similarity between the window centered around it and the windows centered around the other pixels in the whole image, then it takes the similarities as weight to adjust this pixel. It also can be interpreted that the method defines as “neighbor￾hood of a pixel i” any set of pixels j in the image such that a window around j looks like a window around i, and all pixels in that neighborhood are used for pre￾dicting the value at i. In essence, the method uses two windows’ similarity other than two pixels’ to estimate pixels and these windows slide on the whole image. If N2 is the number of pixels in the image, M2 is the size of compare window, the complexity of this al￾gorithm is M2 × N4 . For computational purposes, the simplified algorithm restricts the search of similar win￾dows in a “search window” of size S × S pixels. By fixing the search window at the size of 21 × 21 pixels and the neighborhood size 7×7, the final complexity of the algorithm is 49×441× N2 . However, even the sim￾plified algorithm still takes about 1 minute to denoise a 640 × 480 image on a common PC. It can be seen that the high computational complexity makes it unfeasible to tackle with practical issues. 2.2 Related Work To our knowledge, up to now, there is not much work about nl-means denoising. Lately, Lukin proposed a multiresolution-based approach to improve spatial nl￾means algorithm[16]. The algorithm first performs nl￾means on noisy image in spatial domain with different block sizes and search ranges. Afterward, it transforms the resulting images (generally two images) into fre￾quency domain, and then combines the coefficients by another filter bank, finally inverses the filter bank to get the denoising result. It can be seen that this method only considers getting rid of different frequency noise, but it does not consider the redundancy of image fea￾tures itself and comparing similarity in transform do￾main. Conversely, we compute similarity between im￾age features at different scales, taking advantage of re￾dundancy of different spatial frequency components of image. The methods can not only reduce different fre￾quency noise efficiently but also produce robust results. That is where our approach diverges from Lukin’s al￾gorithm. An algorithm for accelerating spatial nl-means de￾noising was proposed by Mahmoudi and Sapiro[17] . In order to accelerate the algorithm, they introduced filters that eliminate unrelated neighborhoods from weighted average. These filters are based on local aver￾age gray values and gradients, preclassifying neighbor￾hoods and reducing the influence of less-related areas when denoising a given pixel. Note that the times of ac￾celeration of Mahmoudi’s algorithm is related to image content, for example, it is 7 times faster than original spatial nl-means for some images, whereas, for some other images it is more than 20 times faster than the original nl-means. In our method, after speedup, the method is fifty times faster than the original nl-means for all the images. 3 Our Robust and Fast Laplacian Pyramid Based nl-Means 3.1 Laplacian Pyramid Gaussian-Laplacian pyramid has been widely used in texture classification[18,19], feature detection[20] and so on. Given an image I(x, y), the Gaussian pyramid of I is computed as follows: ½ G0(x, y) = I, Gk+1(x, y) = REDUCE(Gk(x, y)), k = 0, 1, . . . , N. (1) The REDUCE operation is carried out by convolving the image with a Gaussian low pass filter and then sub￾sampling the filtered version at a factor of two. Repeat￾ing the REDUCE operation generates Gaussian pyra￾mid. The filter mask is designed such that the center pixel gets more weight than the neighboring ones and the remaining terms are chosen so that their sum is 1. The separable Gaussian kernel is given by: w(r, c) = w(r)w(c). (2) As in the general case, in the paper, the kernel we used is of size 5 × 5, and w(r) = (1/16, 1/4, 3/8, 1/4, 1/16). Thus, the variance of the Gaussian kernel approximates 1.0. In order to generate Laplacian pyramid, the second level of Gaussian image G1 is up-sampled by inserting zeros after each pixel, smoothed with the small kernel as defined in Gaussian pyramid. These interpolation and filtering process is referred to as EXPAND operation. The difference between the previous level of Gaussian image G0 and expanded image Gf1 is named as Lapla￾cian image. This process is continued to obtain a set of band-pass filtered images. Expressed by formulation, Lk(x, y) = Gk(x, y) − EXPAND (Gk+1(x, y)), k = 0, 1, . . . , K, (3)

Yan-Li Liu et al.:A Robust and Fast Non-Local Means Algorithm for Image Denoising 273 where N is the decomposition level and the EXPAND coefficients in the Lk as a weight to estimate (i,) operation is defined as follows. To express as a formula,(i,yi)is: 22 sc.=4m,nc(2,） Lk(红，)= (6) (红：)eLk m=2n=2 k=0,1,,K. (4) where the weight w(i,j)of two coefficients Lk(ri,yi) and Lk(rj,y)depending on their similarity is defined Only terms for which (x-m)/2 and (y-n)/2 are in- as. tegers can be included in the sum. 1 The full reconstruction process is implemented by wk(i,j)=- 间 h (7) first expanding the final level of Gaussian pyramid GK+1,then adding to Lk,k =K,K-1,...,0,thus Here,z(i)is the normalized constant, reconstructing the Gaussian pyramid,level by level,up s(i,) to the original input image,Go.This is a recursive 2(i)= ∑e (8) process,as shown in (5): (jy)ELk Gk Lk +EXPAND(Gk+1),kK,K-1,...,0. Here,hk is a constant proportional to of denoting the (5) noise standard deviation of level k image of Laplacian From the above formulation.it can be seen that pyramid, Laplacian pyramid is a set of bandpass images.It hk ckat, (9) contains most of the image's important textural fea- tures,at different scales.The top level of the pyramid s(i,j)is calculated by the Euclidean distance of the two contains just the highest spatial frequency components coefficients'neighboring region Ni and N;with equal such as the sharp edges,textures,high-frequency noise size (M,M)as: etc.The bottom level contains the lowest spatial fre- quency components.The intermediate levels contain s(亿，)=‖N-N2 (10) features gradually decreasing in spatial frequency from where Ni is a vector formed by concatenating all the high to low. coefficients in the neighborhood of(zi,y). It can be seen that the Laplacian pyramid offers an Since we relate the filter parameters with noise overcomplete representation of the image.Hence,dif- standard deviation of every level of Laplacian pyra- ferent spatial frequency component,e.g.,every level of mid,we need to investigate the noise characteristics Laplacian pyramid is redundant,which makes perform- of Gaussian-Laplacian image pyramid.If the standard ing nl-means on Laplacian image possible.Since the deviation of the original noisy image is o0,the noise features extracted are localized in spatial position and variance o2 of the smoothed image is given by211: scale,they are much more able to express the image, and computing the similarity between them is more rea- 1 sonable and more robust.Meanwhile,Laplacian pyra- G(r;c)drde 4娇o28 (11) mid is translation-invariant,meaning that there is no aliasing effect such as pseudo-Gibbs artifacts in recon- Here,o is the standard deviation of Gaussian kernel. struction. In our paper,o=1.0. Thus,if oo denotes the noise variance of the level 3.2 Our Robust nl-Means 0 image of an Gaussian pyramid,the noise variance of the level k image in Gaussian pyramid is given by Given a noisy image,we first decompose it into a K- 1 level Laplacian pyramid.For convenience,we denote as (o)P=42r6. (12) Lk the level k image of Laplacian pyramid.For a coeffi- cient Lk (i,yi)of Lk (=0,1,...,K),when estimat- The noise variance of every Lk can be computed as fol- ing its denoised version L(i,4),we not only consider lows. the similarities between Lk(i,yi)and its nearby coef- ficients,but also take the similarities between it and (a)2= 1 1 all the coefficients in Lk into account.That is to say, 4r26-4ma26, we exploit the similarity between Lk(i,yi)and all the k=1,2,,N-1. (13)

Yan-Li Liu et al.: A Robust and Fast Non-Local Means Algorithm for Image Denoising 273 where N is the decomposition level and the EXPAND operation is defined as follows. Gk+1(x, y) = 4 X 2 m=2 X 2 n=2 w(m, n)Gk ³x − m 2 , y − n 2 ´ , k = 0, 1, . . . , K. (4) Only terms for which (x − m)/2 and (y − n)/2 are in￾tegers can be included in the sum. The full reconstruction process is implemented by first expanding the final level of Gaussian pyramid GK+1, then adding to Lk, k = K, K − 1, . . . , 0, thus reconstructing the Gaussian pyramid, level by level, up to the original input image, G0. This is a recursive process, as shown in (5): Gk = Lk + EXPAND(Gk+1), k = K, K − 1, . . . , 0. (5) From the above formulation, it can be seen that Laplacian pyramid is a set of bandpass images. It contains most of the image’s important textural fea￾tures, at different scales. The top level of the pyramid contains just the highest spatial frequency components such as the sharp edges, textures, high-frequency noise etc. The bottom level contains the lowest spatial fre￾quency components. The intermediate levels contain features gradually decreasing in spatial frequency from high to low. It can be seen that the Laplacian pyramid offers an overcomplete representation of the image. Hence, dif￾ferent spatial frequency component, e.g., every level of Laplacian pyramid is redundant, which makes perform￾ing nl-means on Laplacian image possible. Since the features extracted are localized in spatial position and scale, they are much more able to express the image, and computing the similarity between them is more rea￾sonable and more robust. Meanwhile, Laplacian pyra￾mid is translation-invariant, meaning that there is no aliasing effect such as pseudo-Gibbs artifacts in recon￾struction. 3.2 Our Robust nl-Means Given a noisy image, we first decompose it into a K￾level Laplacian pyramid. For convenience, we denote as Lk the level k image of Laplacian pyramid. For a coeffi- cient Lk (xi , yi) of Lk (k = 0, 1, . . . , K), when estimat￾ing its denoised version L 0 k (xi , yi), we not only consider the similarities between Lk(xi , yi) and its nearby coef- ficients, but also take the similarities between it and all the coefficients in Lk into account. That is to say, we exploit the similarity between Lk(xi , yi) and all the coefficients in the Lk as a weight to estimate L 0 k (xi , yi). To express as a formula, L 0 k (xi , yi) is: L 0 k (xi , yi) = X (xj ,yj )∈Lk wk(i, j)Lk(xj , yj ) (6) where the weight wk(i, j) of two coefficients Lk(xi , yi) and Lk(xj , yj ) depending on their similarity is defined as, wk(i, j) = 1 z(i) e − s(i,j) h 2 k . (7) Here, z(i) is the normalized constant, z(i) = X (xj ,yj )∈Lk e − s(i,j) h 2 k . (8) Here, hk is a constant proportional to σ k L denoting the noise standard deviation of level k image of Laplacian pyramid, hk = ckσ k L, (9) s(i, j) is calculated by the Euclidean distance of the two coefficients’ neighboring region Ni and Nj with equal size (M, M) as: s(i, j) =k Ni − Nj k 2 (10) where Ni is a vector formed by concatenating all the coefficients in the neighborhood of (xi , yi). Since we relate the filter parameters with noise standard deviation of every level of Laplacian pyra￾mid, we need to investigate the noise characteristics of Gaussian-Laplacian image pyramid. If the standard deviation of the original noisy image is σ 2 0 , the noise variance σ 2 s of the smoothed image is given by[21]: σ 2 s = σ 2 0 Z +∞ −∞ Z +∞ −∞ G 2 σ (r, c)drdc = 1 4πσ2 σ 2 0 . (11) Here, σ is the standard deviation of Gaussian kernel. In our paper, σ=1.0. Thus, if σ 2 0 denotes the noise variance of the level 0 image of an Gaussian pyramid, the noise variance of the level k image in Gaussian pyramid is given by (σ k G) 2 = 1 (4πσ2) k σ 2 0 . (12) The noise variance of every Lk can be computed as fol￾lows. (σ k L) 2 = 1 (4πσ2) k σ 2 0 − 1 (4πσ2) k+1 σ 2 0 , k = 1, 2, . . . , N − 1. (13)

274 J.Comput.Sci.Technol.,Mar.2008,Vol.23,No.2 Obviously,as the level of Laplacian pyramid increases, while N2 and N?can be fast calculated as well us- the amount of noise of Laplacian image decreases ing the SSI scheme proposed in Subsection 3.3.2.Note rapidly. that if the compare window size is Mx M,computing the similarity of the two compare window requires M2 3.3 Acceleration of Our nl-Means Algorithm pixel operations,while,in our algorithm,it is figured out once which is achieved by means of FFT. The most time-consuming part of the nl-means al- gorithm is the calculation of the Euclidian distance be- tween compare windows in image.Our method still need to perform similarity computation for pixels of 防-. every level Lk (=0,1,...,K).In order to speedup the algorithm,we exploit SSI and FFT to transform this process into convolution and summation. 3.3.1 Re-Represent the Similarity of Compare Win- () Image dow For simplicity,we denote L as Lk,k E[0,1,...,N}. Fig.1.Mirrored image. As (10)indicates,the similarity of the pixel L(zi,) and L(xj,y)is computed as: 3.3.2SSl S(i,)=‖N-N2 The principle of SSI resembles Integral image which M-1M-1 has been used in face detection221.For each pixel in = ∑∑L.l,m)-L0m2 the image,integral image maintains the summed value e0m=0 (14) of all the pixels in the upper left part of the original image.Here we extend it to our SSI.Similar to the where Li(l,m)and L;(l,m)represent the correspond- definition of integral image,for each pixel(ro,y0),SSI ing pixels in Ni and Ni respectively. stores its sum for the squared values of the upper left In fact,Li(l,m)in (14)can be represented in the pixels. global coordinates on the mirrored image as:Lj(l- SSI(a0,0）= ∑e (16) j:m-y),with j=3M/2+j,=3M/2+y (see x≤x0,y≤0 Fig.1).So (10)is transformed into: SSI can be obtained in linear time proportional to the image size,we take the following algorithm to cal- M-1M-1 S(i,j)= ∑∑L,m)-L,0-,m-gP culate it efficiently.For zo =0,y0 =0, l=0m=0 SSI(0,0)=L2(0,0): (17) =N?+N-N×N (15) forx0>0,y0=0, where SSI(x0,0)=SSI(x0-1,0)+L2(xo,0): (18) M-1M-1 N2= ∑∑(L,m识， forx0=0,0>0, l=0m=0 M-1M-1 SSI(0,0)=SSI(0,0-1)+L2(0,0): (19) -∑】 2L0-:m-》, forx0>0,y0>0, and SSI(xo,y0)=SSI(zo-1,y0)+SSI(xo,yo-1) M-1M-1 -SSI(x0-1,0-1)+L2(x0,0). N×N=2∑ (L(l,m)·L(l-x,m-) (20) l=0m=0 Obviously,with above algorithm,each pixel in the denotes the convolution between Ni and Ni. original image is processed only once,so the computa- In above formula,Nix Ni can be figured out quickly tional complexity for computing SSI is O(N2),in which with multiplications under the fast Fourier transform, N2 is the size of the image

274 J. Comput. Sci. & Technol., Mar. 2008, Vol.23, No.2 Obviously, as the level of Laplacian pyramid increases, the amount of noise of Laplacian image decreases rapidly. 3.3 Acceleration of Our nl-Means Algorithm The most time-consuming part of the nl-means al￾gorithm is the calculation of the Euclidian distance be￾tween compare windows in image. Our method still need to perform similarity computation for pixels of every level Lk (k = 0, 1, . . . , K). In order to speedup the algorithm, we exploit SSI and FFT to transform this process into convolution and summation. 3.3.1 Re-Represent the Similarity of Compare Win￾dow For simplicity, we denote L as Lk, k ∈ {0, 1, . . . , N}. As (10) indicates, the similarity of the pixel L(xi , yi) and L(xj , yj ) is computed as: S(i, j) = k Ni − Nj k 2 = M X−1 l=0 M X−1 m=0 [Li(l, m) − Lj (l, m)]2 . (14) where Li(l, m) and Lj (l, m) represent the correspond￾ing pixels in Ni and Nj respectively. In fact, Lj (l, m) in (14) can be represented in the global coordinates on the mirrored image as: Lj (l − x 0 j , m − y 0 j ), with x 0 j = 3M/2 + xj , y 0 j = 3M/2 + yj (see Fig.1). So (10) is transformed into: S(i, j) = M X−1 l=0 M X−1 m=0 [Li(l, m) − Lj (l − x 0 j , m − y 0 j )]2 = N 2 i + N 2 j − Ni × Nj (15) where N 2 i = M X−1 l=0 M X−1 m=0 (Li(l, m))2 , N 2 j = M X−1 l=0 M X−1 m=0 (Lj (l − x 0 j , m − y 0 j ))2 , and Ni × Nj = 2 M X−1 l=0 M X−1 m=0 (Li(l, m) · Lj (l − x 0 j , m − y 0 j )) denotes the convolution between Ni and Nj . In above formula, Ni×Nj can be figured out quickly with multiplications under the fast Fourier transform, while N2 i and N2 j can be fast calculated as well us￾ing the SSI scheme proposed in Subsection 3.3.2. Note that if the compare window size is M × M, computing the similarity of the two compare window requires M2 pixel operations, while, in our algorithm, it is figured out once which is achieved by means of FFT. Fig.1. Mirrored image. 3.3.2 SSI The principle of SSI resembles Integral image which has been used in face detection[22]. For each pixel in the image, integral image maintains the summed value of all the pixels in the upper left part of the original image. Here we extend it to our SSI. Similar to the definition of integral image, for each pixel (x0, y0), SSI stores its sum for the squared values of the upper left pixels, SSI(x0, y0) = X x6x0,y6y0 L 2 (x, y). (16) SSI can be obtained in linear time proportional to the image size, we take the following algorithm to cal￾culate it efficiently. For x0 = 0, y0 = 0, SSI(0, 0) = L 2 (0, 0); (17) for x0 > 0, y0 = 0, SSI(x0, 0) = SSI(x0 − 1, 0) + L 2 (x0, 0); (18) for x0 = 0, y0 > 0, SSI(0, y0) = SSI(0, y0 − 1) + L 2 (0, y0); (19) for x0 > 0, y0 > 0, SSI(x0, y0) = SSI(x0 − 1, y0) + SSI(x0, y0 − 1) − SSI(x0 − 1, y0 − 1) + L 2 (x0, y0). (20) Obviously, with above algorithm, each pixel in the original image is processed only once, so the computa￾tional complexity for computing SSI is O(N2 ), in which N2 is the size of the image

Yan-Li Liu et al.:A Robust and Fast Non-Local Means Algorithm for Image Denoising 275 By means of SSI,we can easily get the sum of squares In theory,the neighborhood size M should be congru- for each pixel in any rectangles of the image within con- ent to the size of the repeated patterns in the image. stant time.For example in Fig.2.to calculate the sum In order to achieve convincing results.M should be of squares in rectangle D,only 3 addition operations set larger in higher resolution images,which leads to a are required, significant slowdown in the original algorithm but no performance changing in ours. SD =SAUBUCUD SA-SAUC -SAUB In the paper,we set the decomposition level K=2, =SSI(z1,y1)+SSI(xo,y0)-SSI(xo,1) namely,we decompose every noisy test image into high, -SSI(x1,0). (21) middle and low frequency components.Since there are three images to be processed in our algorithm,at first glance,it seems that our algorithms will be slow.See- ing that the resolution of higher level of pyramid is a quarter of that of the previous level,the window of the B same size in higher level Laplacian pyramid represents larger region of image contents than that of original im- (xo Yo) (6,h) age.Hence,we can restrict the search window in high level of pyramid to a relatively small size.In the paper, we restrict search window in three levels to 21 x 21. (1,y) (xo yo) 11 x 11 and 3 x 3,respectively,from high resolution level to low resolution level.Further,as discussed in Image Section 1,small compare window is sensitive to image details and effective to suppress high frequency noise. Fig.2.Using SSI to compute the summed squared pixels in the Conversely,large compare window is favorable to re- rectangle D. move low frequency noise.Thus we should use small Therefore,N2 and N?can be computed quickly with size compare window on the low level of pyramid and (21)once we get SSI of the noised image.Accurate large size window on the high level of pyramid.In the analysis in following subsection will show that our ac- paper,the compare windows we used are of size 7 x 7, 5x 5.3x 3 respectively,from high resolution level to low celerated algorithm is much faster than the original one. resolution level.Considering resolution factor,compare widows of size 5 x 5 in middle and high level occupy 3.4 Time Complexity Analysis much larger region of original image than that of size 7 x 7 in low level of Laplacian pyramid.Since the size The bottleneck of nl-means is the calculation of the of compare window has no impact on time complexity, Euclidian distance between compare windows in the im- age.In simplified spatial nl-means in which an S x S except the decomposition and synthesis of the Lapla- cian pyramid,the total complexity of our algorithm is search window is used instead of the whole image,it 21×21×N2+11×11×N2/4+3×3×N2/16,ap takes M2x S2 square calculations for each pixel in the proximately 472 x N2.Hence,if we do not count in the image,where M2 is the size of compare window. time of implementing decomposition and synthesis of In our algorithm,this process has been transformed Laplacian pyramid,which is a very fast process due to to computing convolution and summation of squares. its low complexity and simple parallel implementation. It is well known that convolution becomes multiplica our algorithm is about fifty times faster than simplified tion under Fourier transform,thus only S2 multiplica- spatial nl-means. tions need to be taken for each pixel in an image in our algorithm.Furthermore,Fourier transform can be per- formed quickly by modern hardware FFT accelerator 4 Experimental Results and Discussion within a neglectable time,convolution can therefore be fast carried out.As for the summations of squares,only We have tested our method on a set of 8-bit grayscale addition operations are needed in this process which test images.In order to ensure a fair comparison,the can also be neglected compared to the multiplications. test images we used are the same ones used in the de- Thus the computational complexity is Sx S,which is noising experiments reported in [8].We add different M2 times faster than the original algorithm for pixel of standard deviation of additive Gaussian white noise to image.In addition,it can be seen that the change of M the test images,and then make comparisons,in terms has no impact on the time complexity of our method. of both visual quality and PSNR,between the resulting

Yan-Li Liu et al.: A Robust and Fast Non-Local Means Algorithm for Image Denoising 275 By means of SSI, we can easily get the sum of squares for each pixel in any rectangles of the image within con￾stant time. For example in Fig.2, to calculate the sum of squares in rectangle D, only 3 addition operations are required, SD = SA∪B∪C∪D + SA − SA∪C − SA∪B = SSI(x1, y1) + SSI(x0, y0) − SSI(x0, y1) − SSI(x1, y0). (21) Fig.2. Using SSI to compute the summed squared pixels in the rectangle D. Therefore, N2 i and N2 j can be computed quickly with (21) once we get SSI of the noised image. Accurate analysis in following subsection will show that our ac￾celerated algorithm is much faster than the original one. 3.4 Time Complexity Analysis The bottleneck of nl-means is the calculation of the Euclidian distance between compare windows in the im￾age. In simplified spatial nl-means in which an S × S search window is used instead of the whole image, it takes M2 × S 2 square calculations for each pixel in the image, where M2 is the size of compare window. In our algorithm, this process has been transformed to computing convolution and summation of squares. It is well known that convolution becomes multiplica￾tion under Fourier transform, thus only S 2 multiplica￾tions need to be taken for each pixel in an image in our algorithm. Furthermore, Fourier transform can be per￾formed quickly by modern hardware FFT accelerator within a neglectable time, convolution can therefore be fast carried out. As for the summations of squares, only addition operations are needed in this process which can also be neglected compared to the multiplications. Thus the computational complexity is S × S, which is M2 times faster than the original algorithm for pixel of image. In addition, it can be seen that the change of M has no impact on the time complexity of our method. In theory, the neighborhood size M should be congru￾ent to the size of the repeated patterns in the image. In order to achieve convincing results, M should be set larger in higher resolution images, which leads to a significant slowdown in the original algorithm but no performance changing in ours. In the paper, we set the decomposition level K = 2, namely, we decompose every noisy test image into high, middle and low frequency components. Since there are three images to be processed in our algorithm, at first glance, it seems that our algorithms will be slow. See￾ing that the resolution of higher level of pyramid is a quarter of that of the previous level, the window of the same size in higher level Laplacian pyramid represents larger region of image contents than that of original im￾age. Hence, we can restrict the search window in high level of pyramid to a relatively small size. In the paper, we restrict search window in three levels to 21 × 21, 11 × 11 and 3 × 3, respectively, from high resolution level to low resolution level. Further, as discussed in Section 1, small compare window is sensitive to image details and effective to suppress high frequency noise. Conversely, large compare window is favorable to re￾move low frequency noise. Thus we should use small size compare window on the low level of pyramid and large size window on the high level of pyramid. In the paper, the compare windows we used are of size 7 × 7, 5×5, 3×3 respectively, from high resolution level to low resolution level. Considering resolution factor, compare widows of size 5 × 5 in middle and high level occupy much larger region of original image than that of size 7 × 7 in low level of Laplacian pyramid. Since the size of compare window has no impact on time complexity, except the decomposition and synthesis of the Lapla￾cian pyramid, the total complexity of our algorithm is 21 × 21 × N2 + 11 × 11 × N2/4 + 3 × 3 × N2/16, ap￾proximately 472×N2 . Hence, if we do not count in the time of implementing decomposition and synthesis of Laplacian pyramid, which is a very fast process due to its low complexity and simple parallel implementation, our algorithm is about fifty times faster than simplified spatial nl-means. 4 Experimental Results and Discussion We have tested our method on a set of 8-bit grayscale test images. In order to ensure a fair comparison, the test images we used are the same ones used in the de￾noising experiments reported in [8]. We add different standard deviation of additive Gaussian white noise to the test images, and then make comparisons, in terms of both visual quality and PSNR, between the resulting

276 J.Comput.Sci.Technol.,Mar.2008,Vol.23,No.2 images of our algorithm and nl-means algorithm. and the denoised pillar are ragged.In our denoised re- In the paper,as mentioned in Subsection 3.4,we sult,as shown in Fig.4(d),the shape of pillar is well make a 3-level Laplacian pyramid for all the images. preserved.Fig.5 further shows the denoising results Theoretically,we should set different ck for different L with two methods for a natural image.Fig.6 shows the to achieve perfect results.In pursuit of simplicity,we first level and second level of Laplacian images before set the same h(h=10)for all the Lk,and found that and after performing our algorithm for Fig.5(b).Evi- there was no obvious decrease in visual quality between dently,the denoised Laplacian images,see Fig.5(b)and two embodiments. Fig.6(d)perfectly preserve the shape of pillars.Fig.7 is Table 1 summarizes a set of denoising results.As can another comparison of two methods.In the denoising be seen from it,when noise level is not high.our method result of spatial nl-means for "Mandrill"with filtering is slightly better than spatial nl-means in terms of peak parameters h=2.0,as shown in Fig.7(c),the eyes of signal to noise ratio(PSNR).As the noise standard de- mandrill are pale,and on the lower part of image,can- viation increases,the PSNR of our algorithm notably vas effect is pronounced.Fig.7(d)is another result of surpasses that of spatial nl-means.For example,when spatial nl-means with large filtering parameter h=3.5 the noise standard deviation oo is 30,the PSNR of our which is used to remove more noise,it can be seen that method outstrips that of the spatial nl-means with an the resulting image looks more grey than the original excess of 1.0 dB,which can be regarded as a consider- image.In contrast,as shown in Fig.7(e),the eyes of able improvement.In order to better visualize the com- mandrill are quite bright and the image is clean with parison,Fig.3 presents the denoising results as curves, the details well preserved.This is because when noise for the image“Lena” level is very high,in spatial domain,what a compare window represents is not the true content of image ex- Table 1.Denoising Results in PSNR (in [dB]) actly,thus using the similarity of compare window to 60 Lena Peppers estimate pixels would result in large error.In our pa- NL Our Method NL Our Method per,we extract different scale features of image and 10 34.17 34.96 33.32 34.38 operate on these features.The method has two advan- 20 31.45 31.95 31.33 31.87 tages:first,it filters out some noise,thus it is noise 30 28.68 30.08 28.67 29.88 less-sensitive.Second,image features (lines,edges,de- 50 25.41 27.27 26.09 27.36 tails)are more able to represent the image,using it to compute the similarity makes sense and more robust PSNR ■Our Method Spatial Non-Local Means 30 5 00 (a) (b) (c) (d) 01020304050 Fig.3.Comparison of PSNR between our method and spatial Fig.4.Comparison between the cropped result obtained by spa- nl-means. tial nl-means (c)and our method (d)for noisy Lena with oo =20 (b).(a)is the cropped original Lena.(Note that the results of Besides PSNR,we believe that visual quality is also spatial nl-means (c)is smeared out and the denoised polar is an important criterion in judging various noise reduc- ragged.In the result of our method (d),the image is clean and tion methods.Figs.4~7 provide some visual com- the polar is the same as in original image (a).) parisons of results denoised with these two algorithms Fig.4(a)is the cropped part of "Lena"and Fig.4(b) As discussed in Subsection 3.4,our algorithm is is its noisy version with oo =20.We see in Fig.4(c) about fifty times faster than simplified spatial nl-means that the resulting image of spatial nl-means is smeary Experiment on a PC with a Pentium IV 2.4GHz CPU

276 J. Comput. Sci. & Technol., Mar. 2008, Vol.23, No.2 images of our algorithm and nl-means algorithm. In the paper, as mentioned in Subsection 3.4, we make a 3-level Laplacian pyramid for all the images. Theoretically, we should set different ck for different Lk to achieve perfect results. In pursuit of simplicity, we set the same h (h = 10) for all the Lk, and found that there was no obvious decrease in visual quality between two embodiments. Table 1 summarizes a set of denoising results. As can be seen from it, when noise level is not high, our method is slightly better than spatial nl-means in terms of peak signal to noise ratio (PSNR). As the noise standard de￾viation increases, the PSNR of our algorithm notably surpasses that of spatial nl-means. For example, when the noise standard deviation σ0 is 30, the PSNR of our method outstrips that of the spatial nl-means with an excess of 1.0 dB, which can be regarded as a consider￾able improvement. In order to better visualize the com￾parison, Fig.3 presents the denoising results as curves, for the image “Lena”. Table 1. Denoising Results in PSNR (in [dB]) σ0 Lena Peppers NL Our Method NL Our Method 10 34.17 34.96 33.32 34.38 20 31.45 31.95 31.33 31.87 30 28.68 30.08 28.67 29.88 50 25.41 27.27 26.09 27.36 Fig.3. Comparison of PSNR between our method and spatial nl-means. Besides PSNR, we believe that visual quality is also an important criterion in judging various noise reduc￾tion methods. Figs. 4∼7 provide some visual com￾parisons of results denoised with these two algorithms. Fig.4(a) is the cropped part of “Lena” and Fig.4(b) is its noisy version with σ0 = 20. We see in Fig.4(c) that the resulting image of spatial nl-means is smeary and the denoised pillar are ragged. In our denoised re￾sult, as shown in Fig.4(d), the shape of pillar is well preserved. Fig.5 further shows the denoising results with two methods for a natural image. Fig.6 shows the first level and second level of Laplacian images before and after performing our algorithm for Fig.5(b). Evi￾dently, the denoised Laplacian images, see Fig.5(b) and Fig.6(d) perfectly preserve the shape of pillars. Fig.7 is another comparison of two methods. In the denoising result of spatial nl-means for “Mandrill” with filtering parameters h = 2.0, as shown in Fig.7(c), the eyes of mandrill are pale, and on the lower part of image, can￾vas effect is pronounced. Fig.7(d) is another result of spatial nl-means with large filtering parameter h = 3.5 which is used to remove more noise, it can be seen that the resulting image looks more grey than the original image. In contrast, as shown in Fig.7(e), the eyes of mandrill are quite bright and the image is clean with the details well preserved. This is because when noise level is very high, in spatial domain, what a compare window represents is not the true content of image ex￾actly, thus using the similarity of compare window to estimate pixels would result in large error. In our pa￾per, we extract different scale features of image and operate on these features. The method has two advan￾tages: first, it filters out some noise, thus it is noise less-sensitive. Second, image features (lines, edges, de￾tails) are more able to represent the image, using it to compute the similarity makes sense and more robust. Fig.4. Comparison between the cropped result obtained by spa￾tial nl-means (c) and our method (d) for noisy Lena with σ0 = 20 (b). (a) is the cropped original Lena. (Note that the results of spatial nl-means (c) is smeared out and the denoised polar is ragged. In the result of our method (d), the image is clean and the polar is the same as in original image (a).) As discussed in Subsection 3.4, our algorithm is about fifty times faster than simplified spatial nl-means. Experiment on a PC with a Pentium IV 2.4GHz CPU

Yan-Li Liu et al.:A Robust and Fast Non-Local Means Algorithm for Image Denoising 277 (a) (b) (a) (b) (c) (d) (c) (d) Fig.5.Comparison between the results obtained by spatial nl- Fig.6.The second and first level of Laplacian pyramid of Fig.5(b) means and our method for noisy window with oo =30(b).(a) before and after denoising with our method.(a),(c)The second Original image.(b)Noisy window with oo =30.(c)Result by and first level respectively before denoising.(b),(d)Denoised spatial nl-means.(d)Result by our method.It is clearly that second and first level respectively.It can be seen that the noise our method preserves the shape of pillars well. has been removed effectively and the shape of pillars are preserved in (b)and (d) (a) (b) (c) (dΦ (e) Fig.7.Comparison between the results obtained by spatial nl-means and that of our method for noisy Mandrill with oo =30.(a) Original Mandrill image.(b)Noisy Mandrill with o=30.(c)Result obtained by spatial nl-means with filtering parameter h=2.0. (d)Result of the same method with h=3.5.(e)Result by our method.Note that the eyes of Mandrill are pale and the canvas effect appears on the face of Mandrill(c)and (d).In the result of our method (e),the eyes are bright,and there is no artifacts on the face. and 512MB RAM demonstrates that it takes no more requires nearly 10 minutes.This makes the perfor- than 11 seconds to denoise a 500 Megabyte photo using mance of our algorithm acceptable to common users our algorithm,while the original nl-means algorithm as is demonstrated in Table 2

Yan-Li Liu et al.: A Robust and Fast Non-Local Means Algorithm for Image Denoising 277 Fig.5. Comparison between the results obtained by spatial nl￾means and our method for noisy window with σ0 = 30 (b). (a) Original image. (b) Noisy window with σ0 = 30. (c) Result by spatial nl-means. (d) Result by our method. It is clearly that our method preserves the shape of pillars well. Fig.6. The second and first level of Laplacian pyramid of Fig.5(b) before and after denoising with our method. (a), (c) The second and first level respectively before denoising. (b), (d) Denoised second and first level respectively. It can be seen that the noise has been removed effectively and the shape of pillars are preserved in (b) and (d). Fig.7. Comparison between the results obtained by spatial nl-means and that of our method for noisy Mandrill with σ0 = 30. (a) Original Mandrill image. (b) Noisy Mandrill with σ0 = 30. (c) Result obtained by spatial nl-means with filtering parameter h = 2.0. (d) Result of the same method with h = 3.5. (e) Result by our method. Note that the eyes of Mandrill are pale and the canvas effect appears on the face of Mandrill (c) and (d). In the result of our method (e), the eyes are bright, and there is no artifacts on the face. and 512MB RAM demonstrates that it takes no more than 11 seconds to denoise a 500 Megabyte photo using our algorithm, while the original nl-means algorithm requires nearly 10 minutes. This makes the perfor￾mance of our algorithm acceptable to common users as is demonstrated in Table 2

278 J.Comput.Sci.Technol.,Mar.2008,Vol.23,No.2 (5]Donoho D.De-noising by soft-thresholding.IEEE Trans.In- Table 2.Performance Results formation Theory,1995,41(3):613-627. [6]Chambolle A,DeVore R A,Lee N Y,Lucier B J.Nonlinear Image Size Simplified Our Fast Ratio wavelet image processing:Variational problems,compression, Spatial NL (s)Method (s) and noise removal through wavelet shrinkage.IEEE Trans. 512×512 28.16 0.38 74.1 Image Processing,1998,7(1):319-335. 1024×768 85.45 1.54 55.5 [7]Cohen I,Raz S,Malah D.Translation invariant denoising 2592×1944 548.1 10.63 51.6 using the minimum description length criterion.Signal Pro- cessing,1999,75(3):201-223. [8 Portilla J,Strela V,Wainwright M J,Simoncelli E P.Image denoising using scale mixtures of Gaussians in the wavelet 5 Conclusion domain.TEEE Trans.Image Processing,2003,12(11):1338- 1351. In the paper,we proposed a robust and fast non- [9]Romberg J,Choi H,Baraniuk R G.Bayesian tree-structured local denoising algorithm.The algorithm is based on wavelet-domain image modeling using hidden Markov models. Laplacian pyramid.Due to the redundancy property IEEE Trans.Image Processing,2001,10(7):1056-1068. of Laplacian pyramid,we can further perform nl-means [10]Buades A,Coll B,Morel J M.A non-local algorithm for image building upon image redundancy in Laplacian pyramid. denoising.In Proc.IEEE Computer Society Conference on Computer Vision and Pattern Recognition.San Diego,USA. We employ Laplacian pyramid to break up noisy im- 2005,Vol.2.pp.60-65. age into bandpassed images.By performing nl-means [11]Ahmadian A,Bharath AA.Orthogonal wavelets for image on different levels of Laplacian pyramid with differ- transmission and compression schemes:Implementation and ent sizes of compare windows,we effectively remove results.In Proc.SPIE,1996,2825(2):822-833. high-frequency noise and low-frequency noise while pre- [12]Kharate G K,Ghatol AA,Rege PP.Image compression us- serving the image details (edges,textures,etc.).We ing wavelet packet tree based on threshold entropy.In Proc. the 24th IASTED International Conference on Signal Pro- also utilize SSI and FFT to propose a speedup algo- cessing,Pattern Recognition,and Applications,Innsbruck, rithm.After acceleration,our method is nearly fifty Austria,2006,pp.322-325. times faster that spatial nl-means algorithm.It can be [13 Jansen M,Bultheel A.Multiple wavelet threshold estima- seen that our accelerated algorithm is more feasible to tion by generalized crossvalidation for images with correlated noise.IEEE Trans.Image Processing,1999,8(7):947-953. tackle with practical problems. [14]Aiazzi B,Alparone L,Baronti S,Borri G.Pyramid-based mul- Compared with spatial nl-means results,the result- tiresolution adaptive filters for additive multiplicative image ing image of our algorithm is neat,while spatial nl- noise.IEEE Trans.Circuits Syst.II,1998,45(8):1092- means results are rather blotchy.Our result is a bit 1097. more blurry but introduces less artifacts such as can- [15 Burt P J,Adelson E H.The Laplacian pyramid as a com- vas effect than spatial nl-means.The limitation of the pact image code.IEEE Trans.Communications,1983,31(4): 532-540. method,as shown in Table 2,is:when images to be [16]Alexey L A.Multiresolution approach for improving qual- denoised reach to a particularly large size,the ratio be ity of image denoising algorithms.In Proc.IEEE Interna- tween times needed in our method and spatial nl-means tional Conference on Acoustics,Speech and Signal Process- algorithms slowly decrease.This is because,for large ing,Toulouse,France,2006. image,it takes several minutes to decompose it into [17]Mahmoudi M,Sapiro G.Fast image and video denoising via Laplacian pyramid.We suggest this limitation can be nonlocal means of similar neighborhoods.TEEE Signal Pro- cessing Letters,2005,12(12):839-842. further overcomed by using GPU to implement the de- [18 Heeger D J,Bergen J R.Pyramid-based texture analy- composition process. sis/synthesis.In Proc.SIGGRAPH,Los Angeles,USA,1995. pp.229-238. References [19]Greenspan H,Goodman R,Chellappa R,Anderson C H. Learning texture discrimination rules in a multiresolution sys- tem.IEEE Trans.Pattern Analysis and Machine Intelli- [1]Lindenbaum M,Fischer M,Bruckstein A M.On Gabor con- 9 gence,.1994,16(9):894-901. tribution to image enhancement.Pattern Recognition,1994, [20]Bouzidi A,Baaziz N.Contourlet domain feature extraction for 27(1):1-8. [2]Alvarez L,Lions P L,Morel J M.Image selective smooth- image content authentication.In Proc.IEEE 8th Workshop ing and edge detection by nonlinear diffusion (ii).Journal of on Multimedia Signal Processing,Victoria,Canada,2006, Numerical Analysis,1992,29(3):845-866. pp.202-206. [3]Yin L,Yang R,Gabbouj M,Neuvo Y.Weighted median fil- 21]Fuchs C.Extraktion polymorpher Bildstrukturen und ihre ters:A tutorial.IEEE Trans.Circuits and Systems,1996, topologische und geometrische Gruppierung.DGK,Bayer. 43(3):157-192. Akademie der Wissenschaften,Reihe C,Heft 502,1998. [4]Tomasi C,Manduchi R.Bilateral filtering for gray and color [22]Viola P,Michael J.Rapid object detection using a boosted images.In Proc.the Sirth International Conference on Com- cascade of simple features.In Proc.IEEE CVPR,Hanoii, puter Vision,Bombay,India,1998,pp.839-846. USA,2001,pp.511-518

278 J. Comput. Sci. & Technol., Mar. 2008, Vol.23, No.2 Table 2. Performance Results Image Size Simplified Our Fast Ratio Spatial NL (s) Method (s) 512 × 512 28.16 0.38 74.1 1024 × 768 85.45 1.54 55.5 2592 × 1944 548.1 10.63 51.6 5 Conclusion In the paper, we proposed a robust and fast non￾local denoising algorithm. The algorithm is based on Laplacian pyramid. Due to the redundancy property of Laplacian pyramid, we can further perform nl-means building upon image redundancy in Laplacian pyramid. We employ Laplacian pyramid to break up noisy im￾age into bandpassed images. By performing nl-means on different levels of Laplacian pyramid with differ￾ent sizes of compare windows, we effectively remove high-frequency noise and low-frequency noise while pre￾serving the image details (edges, textures, etc.). We also utilize SSI and FFT to propose a speedup algo￾rithm. After acceleration, our method is nearly fifty times faster that spatial nl-means algorithm. It can be seen that our accelerated algorithm is more feasible to tackle with practical problems. Compared with spatial nl-means results, the result￾ing image of our algorithm is neat, while spatial nl￾means results are rather blotchy. Our result is a bit more blurry but introduces less artifacts such as can￾vas effect than spatial nl-means. The limitation of the method, as shown in Table 2, is: when images to be denoised reach to a particularly large size, the ratio be￾tween times needed in our method and spatial nl-means algorithms slowly decrease. This is because, for large image, it takes several minutes to decompose it into Laplacian pyramid. We suggest this limitation can be further overcomed by using GPU to implement the de￾composition process. References [1] Lindenbaum M, Fischer M, Bruckstein A M. On Gabor con￾tribution to image enhancement. Pattern Recognition, 1994, 27(1): 1–8. [2] Alvarez L, Lions P L, Morel J M. Image selective smooth￾ing and edge detection by nonlinear diffusion (ii). Journal of Numerical Analysis, 1992, 29(3): 845–866. [3] Yin L, Yang R, Gabbouj M, Neuvo Y. Weighted median fil￾ters: A tutorial. IEEE Trans. Circuits and Systems, 1996, 43(3): 157–192. [4] Tomasi C, Manduchi R. Bilateral filtering for gray and color images. In Proc. the Sixth International Conference on Com￾puter Vision, Bombay, India, 1998, pp.839–846. [5] Donoho D. De-noising by soft-thresholding. IEEE Trans. In￾formation Theory, 1995, 41(3): 613–627. [6] Chambolle A, DeVore R A, Lee N Y, Lucier B J. Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage. IEEE Trans. Image Processing, 1998, 7(1): 319–335. [7] Cohen I, Raz S, Malah D. Translation invariant denoising using the minimum description length criterion. Signal Pro￾cessing, 1999, 75(3): 201–223. [8] Portilla J, Strela V, Wainwright M J, Simoncelli E P. Image denoising using scale mixtures of Gaussians in the wavelet domain. IEEE Trans. Image Processing, 2003, 12(11): 1338– 1351. [9] Romberg J, Choi H, Baraniuk R G. Bayesian tree-structured wavelet-domain image modeling using hidden Markov models. IEEE Trans. Image Processing, 2001, 10(7): 1056–1068. [10] Buades A, Coll B, Morel J M. A non-local algorithm for image denoising. In Proc. IEEE Computer Society Conference on Computer Vision and Pattern Recognition, San Diego, USA, 2005, Vol.2, pp.60–65. [11] Ahmadian A, Bharath A A. Orthogonal wavelets for image transmission and compression schemes: Implementation and results. In Proc. SPIE, 1996, 2825(2): 822–833. [12] Kharate G K, Ghatol A A, Rege P P. Image compression us￾ing wavelet packet tree based on threshold entropy. In Proc. the 24th IASTED International Conference on Signal Pro￾cessing, Pattern Recognition, and Applications, Innsbruck, Austria, 2006, pp.322–325. [13] Jansen M, Bultheel A. Multiple wavelet threshold estima￾tion by generalized crossvalidation for images with correlated noise. IEEE Trans. Image Processing, 1999, 8(7): 947–953. [14] Aiazzi B, Alparone L, Baronti S, Borri G. Pyramid-based mul￾tiresolution adaptive filters for additive multiplicative image noise. IEEE Trans. Circuits Syst. II, 1998, 45(8): 1092– 1097. [15] Burt P J, Adelson E H. The Laplacian pyramid as a com￾pact image code. IEEE Trans. Communications, 1983, 31(4): 532–540. [16] Alexey L A. Multiresolution approach for improving qual￾ity of image denoising algorithms. In Proc. IEEE Interna￾tional Conference on Acoustics, Speech and Signal Process￾ing, Toulouse, France, 2006. [17] Mahmoudi M, Sapiro G. Fast image and video denoising via nonlocal means of similar neighborhoods. IEEE Signal Pro￾cessing Letters, 2005, 12(12): 839–842. [18] Heeger D J, Bergen J R. Pyramid-based texture analy￾sis/synthesis. In Proc. SIGGRAPH, Los Angeles, USA, 1995, pp.229–238. [19] Greenspan H, Goodman R, Chellappa R, Anderson C H. Learning texture discrimination rules in a multiresolution sys￾tem. IEEE Trans. Pattern Analysis and Machine Intelli￾gence, 1994, 16(9): 894–901. [20] Bouzidi A, Baaziz N. Contourlet domain feature extraction for image content authentication. In Proc. IEEE 8th Workshop on Multimedia Signal Processing, Victoria, Canada, 2006, pp.202–206. [21] Fuchs C. Extraktion polymorpher Bildstrukturen und ihre topologische und geometrische Gruppierung. DGK, Bayer. Akademie der Wissenschaften, Reihe C, Heft 502, 1998. [22] Viola P, Michael J. Rapid object detection using a boosted cascade of simple features. In Proc. IEEE CVPR, Hanoii, USA, 2001, pp.511–518

Yan-Li Liu et al.:A Robust and Fast Non-Local Means Algorithm for Image Denoising 279 Yan-Li Liu received her B.S.de- Yan-Wen Guo received his gree in applied mathematics from Ph.D.degree in applied mathemat- Henan University in 2004.Now,she ics from State Key Lab of CAD&CG, is a Ph.D.candidate in State Key Zhejiang University in 2006.He is Lab of CAD&CG,Zhejiang Univer- currently an assistant professor at sity,majoring in applied mathemat- the National Laboratory of Novel ics.Her research interests include Software Technology,Nanjing Uni- image and video processing,face re- versity.He worked as a visiting re- lated research and augmented reality. searcher in Department of Computer Science and Engineering,the Chinese University of Hong Kong in 2006.His main research inter- Jin Wang received his Ph.D.de ests include image and video processing,geometry process- gree in computer science and technol- ing,and face related applications. ogy from Zhejiang University,P.R. China in 2007.Currently he is Qun-Sheng Peng graduated a researcher in State Key Lab of from Beijing Mechanical College in CAD&CG.His research interests fo- 1970 and received a Ph.D.degree cus on photo/video enhancement, from the Department of Comput- face relevant research,computer vi- ing Studies,University of East An- sion&computer graphics. glia,UK in 1983.Now he is a pro- fessor of computer graphics at Zhe- jiang University.He serves currently Xi Chen received her B.S.degree as a member of the editorial boards in computer science and technology of several international and Chinese from Zhejiang University,P.R.China journals.His research interests include realistic image syn- in 2007.She is currently pursuing her thesis,computer animation,scientific data visualization, M.Sc.degree in Department of Com virtual reality,bio-molecule modeling.He is a senior mem- puter Science,University of British ber of CCF. Columbia.Canada.Her research in- terests include image and video pro- cessing

Yan-Li Liu et al.: A Robust and Fast Non-Local Means Algorithm for Image Denoising 279 Yan-Li Liu received her B.S. de￾gree in applied mathematics from Henan University in 2004. Now, she is a Ph.D. candidate in State Key Lab of CAD&CG, Zhejiang Univer￾sity, majoring in applied mathemat￾ics. Her research interests include image and video processing, face re￾lated research and augmented reality. Jin Wang received his Ph.D. de￾gree in computer science and technol￾ogy from Zhejiang University, P.R. China in 2007. Currently he is a researcher in State Key Lab of CAD&CG. His research interests fo￾cus on photo/video enhancement, face relevant research, computer vi￾sion& computer graphics. Xi Chen received her B.S. degree in computer science and technology from Zhejiang University, P.R. China in 2007. She is currently pursuing her M.Sc. degree in Department of Com￾puter Science, University of British Columbia, Canada. Her research in￾terests include image and video pro￾cessing. Yan-Wen Guo received his Ph.D. degree in applied mathemat￾ics from State Key Lab of CAD&CG, Zhejiang University in 2006. He is currently an assistant professor at the National Laboratory of Novel Software Technology, Nanjing Uni￾versity. He worked as a visiting re￾searcher in Department of Computer Science and Engineering, the Chinese University of Hong Kong in 2006. His main research inter￾ests include image and video processing, geometry process￾ing, and face related applications. Qun-Sheng Peng graduated from Beijing Mechanical College in 1970 and received a Ph.D. degree from the Department of Comput￾ing Studies, University of East An￾glia, UK in 1983. Now he is a pro￾fessor of computer graphics at Zhe￾jiang University. He serves currently as a member of the editorial boards of several international and Chinese journals. His research interests include realistic image syn￾thesis, computer animation, scientific data visualization, virtual reality, bio-molecule modeling. He is a senior mem￾ber of CCF