基于小波預處理及自適應權函數的非局部均值濾波算法

毛玉星,李 超,賈海威

(重慶大學 電氣工程學院,重慶 400030)

1 引 言

圖像去噪技術一直以來都是圖像處理領域的熱點與關鍵技術,目的是從含噪圖像得到真實圖像信息,退化模型如(1)所示,

Y=X+n

(1)

其中,Y是含噪圖像矩陣,也可稱為觀測圖像,X是真實圖像矩陣,n為加性噪聲,在本文假設n為加性高斯白噪聲,即n～N(0,σ2).圖像去噪的目的是由觀測圖像Y通過變換計算,最大程度的準確獲得X.為了準確的獲取真實圖像,主要研究的空域濾波算法有均值濾波、高斯濾波、非局部均值濾波[1,2]、雙邊濾波[3]等,頻域濾波算法有維納濾波、傅里葉變換濾波、小波閾值濾波等[4,5].其中,非局部均值濾波(NLM)因良好的去噪效果受到了廣泛關注,其充分利用圖像的冗余信息,以欲消噪點為中心,在一定搜索范圍內,以每一像素點與中心點的加權歐式距離大小為權重,對周圍點加權求和來得到中心點.

非局部均值濾波算法的相關研究主要集中在兩個方面:1)對算法速度的提升,因其在計算歐氏距離過程中,是點與點的計算與平移,存在大量的計算,時間較長.為了優化加速,引入積分圖的思想,先將含噪圖像偏移,并與原圖做差分平方,作為圖像的歐氏距離矩陣,再對其求積分得到積分矩陣[6-8].此時,計算任意兩個鄰域框的歐氏距離時,只需要以待復原點為中心的鄰域框四個角點在積分矩陣中的值即可,大幅度縮短了計算時間,為實時處理提供了方法,但是該方法使用圖像塊的思想,不能根據圖像的局部信息適時調整濾波參數,對于圖像的去噪效果有一定影響;2)對算法去噪效果的提升,主要是中心點權重、鄰域框尺寸、濾波平滑參數等的優化及自適應.對于中心像素的權重,在文獻[9-11]中,作者分別使用斯坦無偏風險估計(SURE)、最大似然估計、James-Stein 中心收縮估計來確定中心點權重,相較于通常讓中心點權重等于搜索范圍內的最大權重或者直接令中心點權重為1,在去噪效果上有一定程度的改善,但是會增加運算的復雜度.對于鄰域框的研究,主要是對圖像邊緣處鄰域框尺寸的選擇和非相關點的濾除,在文獻[12,13]中,作者使用結構張量矩陣來表示圖像紋理結構的復雜度,將其分級,不同復雜等級的像素點,使用不同大小的鄰域框,在文獻[14,15]中,根據不同的理論建立權重閾值,當搜索范圍內匹配點的權重小于閾值的時候,則認為其為不相關點,使其權重為0.對于濾波平滑參數,一般認為其與圖像的噪聲方差成正比,但同時也與圖像局部結構信息相關,在文獻[16]中,作者使用SUSAN邊緣檢測理論,以邊緣檢測響應來修正原始平滑參數(通過噪聲方差得到),提升濾波效果,在文獻[17]中,由均方誤差來評價去噪效果,建立關于平滑參數的函數,再通過黃金分割算法計算最佳平滑參數以使均方誤差最小,達到自適應最佳平滑參數的目的.對于歐氏距離權函數,一般使用高斯權函數,有的研究中也使用矩形權函數[18],但現階段對于高斯權函數和矩形權函數在非局部均值濾波中的區別和對去噪結果的影響,研究很少.此外,也有通過對圖像進行預處理提升去噪效果的研究,在文獻[19]中,作者提出使用兩步去噪的方法,第一步用較小的平滑參數hbasic進行非局部均值濾波作為對圖像的預處理,第二步用最佳平滑參數hfinal進行非局部均值濾波,能夠較準確的得到去噪圖像,雖然該方法能夠提高去噪效果,但是耗時長,且會在第一步處理時對圖像邊緣產生平滑.

通過上述分析,本文提出使用小波閾值去噪作為非局部均值濾波的預處理,然后用自適應歐氏距離權函數的非局部均值濾波進行去噪.后續章節分布如下:第二部分對非局部均值濾波算法作了簡單回顧,第三部分引入適用于預處理的小波閾值去噪方法,第四部分探究了歐氏距離權函數與圖像噪聲強度以及局部結構的關系,提出自適應權函數選擇方法,第五部分對本文所提方法進行實驗仿真,第六部分對全文進行總結.

2 非局部均值濾波算法

含噪圖像Y中任意像素點y(i,j),運用非局部均值濾波算法去噪得x(i,j),如公式(2)所示[6]:

(2)

其中,Ω為以y(i,j)為中心的圖像二維搜索空間,y(p,q)為該空間內與y(i,j)不重合的任意點.W(i,j,p,q)為y(p,q)的權重,計算公式如公式(3)所示,w(i,j)為y(i,j)的權重,一般選用搜索框范圍內的最大權重,z為歸一化參數,計算如公式(4)所示.

(3)

(4)

其中,M為一個二維的權函數,一般使用高斯權函數或者矩形權函數,*為卷積運算,‖·‖2為L2范數,P(i,j)和P(p,q)是指以y(i,j)和y(p,q)為中心的鄰域塊,也叫相似框.Di,j為相似框歐式距離矩陣與M的卷積,即加權歐式距離,h為濾波平滑參數.

3 小波閾值去噪預處理

小波分析作為信號處理的工具,是繼Fourier分析之后又一有效的時頻分析方法,可同時進行時域和頻域分析,且有時域局部化和多分辨率特性,最早由Weaver將小波變換用于圖像降噪[4].所謂的小波閾值降噪就是先將圖像進行小波分解,再將小波系數的幅值同一個閾值進行比較,若小波系數的幅值比這個閾值小,則把小波系數置為0,若這個幅值比小波閾值大,則把小波系數進行保留和修改,最后將小波系數進行重構得到去噪圖像[20,21].

小波閾值去噪分為硬閾值法和軟閾值法,分別如公式(5)、(6)所示.

(5)

(6)

式中,T為閾值,y為含噪圖像的小波變換系數,sgn(y)表示y的符號,Thard、Tsoft分別為硬閾值法濾波和軟閾值法濾波的收縮函數.因為軟閾值法中小波系數的幅值被減去了一部分,因此軟閾值濾波中小波參數的估計是有偏的,濾波后的信號會過于平滑,丟失圖像細節信息.因此,在本文中使用硬閾值法進行濾波,使用統一閾值Tuniv,如公式(7)所示.

(7)

式中,σn為零均值加性高斯白噪聲的標準差,可由小波變換系數估計,如公式(8)所示,N為小波系數的總個數.

(8)

本文中,使用小波閾值去噪的目的是為非局部均值濾波進行預處理,即在頻域濾除一部分噪聲,來提高NLM的效果,所以在小波去噪時要盡可能的保留圖像的結構以及邊緣信息.通過測試,使用小波基函數sym8在預處理中有最好的表現.此處使用分辨率為256*256,噪聲標準差σ=25的8位cameraman灰度圖,對圖像單層小波分解與雙層分解的閾值去噪結果進行對比,如圖1所示.由圖可知,雙層分解去噪能夠獲得更好的去噪效果,但是圖像更加模糊,濾除了大量的結構信息,這必然會對后續步驟中的NLM帶來不利影響,因此,本文使用單層分解硬閾值小波去噪進行預處理.

圖1 小波閾值去噪結果Fig.1 Wavelet threshold denoising

4 自適應權函數的研究

在進行兩個鄰域塊P(i,j)和P(p,q)的加權歐式距離計算時,一般使用高斯權函數模板與所得歐氏距離矩陣卷積求和[1,2].然而,研究發現,在高噪聲的圖像和低噪聲的圖像平滑區域使用矩形權函數模板能夠取得更好的去噪效果,因此,有的研究中直接使用矩形權函數[6],但是對于權函數的自適應選擇研究較少.權函數模板與相似框有一樣的尺寸,高斯權函數和矩形權函數的計算公式分別見公式(9)與(10).

(9)

K(i,j)=1/(2p+1)2

(10)

式中,α為高斯權函數標準差,p為相似框半徑.由公式(9)可知,權函數標準差越小,計算所得高斯模板中心點的值越大,即在加權歐氏距離計算中,中心點所占權重越大;而當權函數標準差越大,高斯模板中各數值相差越小,直到所有值趨于相等,變成公式(10)計算的矩形模板.在傳統NLM中對高斯和矩形模板的選擇,以及高斯權函數標準差的選擇一般由經驗隨機選擇,且一旦選擇之后,對整幅圖的每一個像素點都是采用相同的權函數.但是在去噪過程中,因圖像紋理結構以及含噪強度的不同,相似框歐式距離計算時中心點的權重占比應該有所區別.本文以灰度圖為例分析,對無噪圖像,兩個像素點相似性判斷只計算這兩個點的歐式距離即可,即相似框的中心點權重為1,周圍點權重全為0,而當圖像噪聲強度越大,則越依賴相似框中周圍點的歐氏距離來輔助判斷,中心點的權重越小.因此,在噪聲強度低時,應使用標準差較小的高斯模板,隨著噪聲增加,高斯權函數標準差增大,以加大相似框中非中心點的權重,直到權函數模板趨于矩形模板,此時,相似框中所有像素點的權重相等.對圖像局部而言,在低頻平坦區,相似框中每個點灰度值相差不大,所有點的權重應該相近,使用矩形模板較好,而在高頻邊緣區,中心點與周圍點的灰度值相差較大,應該減小周圍點權重,使用標準差較小的高斯權函數模板.

4.1 圖像局部分析

由上述分析可知,權函數的使用與圖像局部紋理結構有關.為了更加準確的描述圖像的局部紋理復雜度,本文使用線性結構張量對其進行分類[12],線性結構張量由初始結構張量與高斯核函數卷積得到,如公式(11)所示.

(11)

其中Gδ為標準差為δ的高斯核函數,Ix、Iy分別為圖像水平方向與垂直方向的梯度.對公式(11)進行矩陣變換,可求得法線方向特征值λ1和切線方向特征值λ2如公式(12)所示.

(12)

定義λ=|λ1(i,j)-λ2(i,j)|,λ(i,j)越大,則該點灰度變化越劇烈,反之,則圖像該點處越平坦.可據此將圖像的紋理復雜度分為4級,并根據不同的等級使用不同的權函數M,如公式(13)所示,K(f)表示半徑為f的矩形模板,G(α,f)表示標準差為α,半徑為f的高斯模板.

(13)

其中,α1>α2>α3,n1

4.2 圖像整體分析

由上述分析可知,權函數的選擇與圖像噪聲強度相關.通過實驗發現,當使用高斯權函數且當α>3時,與使用矩形權函數去噪的峰值信噪比(PSNR)相近,即高斯權函數模板近似為矩形權函數模板.本文分別對預處理之后的Lena圖、Peppers圖、boat圖和cameraman圖進行實驗,探究圖像噪聲強度與權函數的關系,以及在不同噪聲情況下的最佳高斯權函數標準差α.使用高斯權函數進行測試,權函數標準差間隔0.2,從0.6一直到3,取不同噪聲強度時最佳權函數標準差如圖2所示,從圖中可以看出,在圖像噪聲σ≥35時高斯權函數最佳標準差穩定在最大值3,此時高斯模板近似為矩形模板.然后,為了探究高斯權函數與矩形權函數的適用范圍,將不同噪聲情況下使用高斯權函數去噪的最佳PSNR和使用矩形權函數去噪的PSNR做差如圖3所示,由圖可知,圖像含噪越少,使用高斯權函數效果越優于矩形權函數,在圖像噪聲標準差為10時,cameraman圖使用高斯權函數去噪的PSNR比使用矩形權函數高0.8dB.隨著噪聲強度變大,曲線趨于負,結合圖2分析,之所以負值時幅值較小是因為高斯權函數標準差為3,此時的高斯模板已近似為矩形模板,故高噪聲時使用矩形權函數效果更優.綜上,當σ≥35時,使用矩形權函數,σ<35時,使用高斯權函數,且標準差設置為α=(0.08～0.1)×σ,α1=α,α2=0.8×α,α3=0.6×α往往能夠取得最好的去噪效果,此處取n1為0.05,n2為0.3,n3為0.6.

圖2 權函數最佳標準差與圖像噪聲關系Fig.2 Relationship between the kernel function′s optimal standard deviation and the image′s noise level

圖3 兩種權函數PSNR差值與圖像噪聲關系Fig.3 Relationship between the difference of PSNR calculated by two different kernel function and the image′s noise level

綜上分析可知,本文對一幅圖像使用非局部均值去噪步驟如下:1)估計圖像噪聲方差,可由公式(8)估計,本文假設圖像噪聲方差已知;2)對圖像進行小波單層硬閾值去噪預處理;3)根據噪聲強度自適應選擇權函數對預處理之后的圖像非局部均值濾波,權函數選擇如公式(14)所示,G表示高斯權函數,K表示矩形權函數,當選擇G,再結合公式(13)根據局部結構自適應設置高斯權函數標準差.

(14)

5 仿真與實驗分析

本文在matlab環境中使用C語言編程進行仿真,選取搜索框尺寸S=21×21,相似框尺寸P=5×5,濾波參數h從0.6*σ到σ間距0.5取每次實驗最佳的PSNR,原始NLM中高斯權函數標準差隨機選取1.4.為了驗證本文中小波預處理以及自適應權函數各自的作用,以Lena圖和Peppers圖為例,在不同噪聲強度下,對比文獻[1]中原始NLM、小波預處理之后分別使用高斯權函數和矩形權函數以及本文自適應權函數進行非局部均值濾波的PSNR見圖4.其中,2-G、2-K分別表示小波預處理之后使用高斯權函數和矩形權函數進行非局部均值濾波,高斯權函數的標準差選擇1.4.對比2-G和文獻[1],可以看出小波預處理對去噪結果的PSNR提高很大,在σ=35時,對Peppers圖高達0.7dB,且隨著噪聲強度變大,有進一步提高的趨勢;對比2-G和2-K曲線可以發現,噪聲低時使用高斯權函數效果優于使用矩形權函數,而當噪聲強度大時,使用矩形權函數有更好的去噪效果,驗證了根據不同噪聲強度選用不同的權函數的必要.本文的自適應權函數結合了2-G和2-K兩種權函數的優點,而本文方法在低噪聲時效果優于2-G是因為本文方法根據噪聲強度和圖像紋理結構,自適應確定高斯權函數的標準差.

圖4 使用不同去噪方法的峰值信噪比(左邊Lena圖、右邊Peppers圖)Fig.4 PSNR comparison with different method(Left is Lena,right is Peppers)

表1 不同方法去噪的峰值信噪比(分貝)
Table 1 PSNR results of various methods(dB)

σ1520253035405060Lena 256?256噪聲圖24.692922.293420.374418.90417.651616.569914.82413.4823文獻131.372729.771628.418827.370726.344125.432824.002422.7428文獻1931.098730.036528.968527.468926.532725.645524.499823.2198本文方法31.581630.170929.073528.090627.225726.319425.185124.0877Peppers 256?256噪聲圖24.664822.261220.33118.752117.479516.459314.689813.3302文獻131.857430.2628.872727.533326.530225.609924.186622.8341文獻1931.572730.54929.571927.881626.946925.980224.779123.4415本文方法32.185230.851229.647428.456827.557226.738325.538524.2277

圖5 不同方法的消噪圖對比,(a)、(b)、(c)、(d)、(e)分別是原圖、噪聲圖(Lena、Peppers噪聲標準差分別為20、25)、原始NLM方法去噪圖、文獻[19]兩步去噪圖、本文方法去噪圖Fig.5 Filtering performance of various denoising method(a)Oringal Image(b)Noisy Image,Lena′s σ=20,Peppers′s σ=25(c)Oringal NLM(d)Two Stage Method[19](e)Our Method

將兩種圖像運用不同方法去噪的PSNR統計如表1所示.從表中數據分析可知,相較于前兩種方法,本文方法有更高的PSNR,在噪聲強度越大時,提升效果越明顯,與原始NLM相比,PSNR最高有1.4dB的提升.將Lena圖 =20,Peppers圖=25時運用不同方法去噪的效果圖展示如圖5,其中圖5(a)為原圖、圖5(b)為噪聲圖、圖5(c)為原始NLM效果圖、圖5(d)為文獻[19]方法的效果圖、圖5(e)為本文方法消噪效果圖.從圖5可以看出,運用本文方法進行去噪,在圖像邊緣輪廓等細節處有更好的視覺效果.原因是對相同噪聲強度的圖像,經過小波預處理,在進行非局部均值濾波時相比于原始NLM選擇的最佳平滑參數更小,即圖像平滑程度更低,且本文根據圖像邊緣結構自適應選擇權函數,能更好的保留邊緣結構.

表2 不同去噪方法時間(秒)
Table 2 Time consumption of different method(s)

NLM文獻[19]本文方法256?25663.93s127.66s64.93s512?512253.25s509.77s258.9s

將不同分辨率的Lena含噪圖運用不同去噪方法的時間對比如表2所示,從表中可以看出,本文相比于文獻[19]中兩步濾波方法,同樣進行了預處理,但是時間節省了近一半,相比于原始算法,本文方法在時間增加很少的情況下,提高了去噪的效果.

6 總 結

通過對非局部均值濾波和小波閾值去噪的分析,本文提出將非局部均值濾波與小波去噪有效結合的方法,對一幅含噪圖像首先通過小波閾值去噪預處理,再使用非局部均值濾波.這種方法可以在不需要優化非局部均值濾波參數的情況下,一定程度的提高去噪效果,圖像噪聲強度越大時提升越明顯,且相對于其他預處理方式,使用小波預處理可以更好的保存圖像邊緣等結構信息.此外,本文對非局部均值濾波的歐氏距離權函數進行了研究,發現在低噪聲圖像的高頻邊緣區使用高斯權函數,高噪聲圖像和低噪聲圖像的平坦區使用矩形權函數相較于原始算法有更好的去噪效果,提出基于圖像噪聲強度與紋理結構的自適應權函數非局部均值濾波.最后,經過仿真對比驗證,本文方法有更高的峰值信噪比和更好的視覺效果.