计算机科学与技术（参考文献）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