改進非局部均值各向異性擴散圖像去噪算法*

王 磊，王 敏，張鵬程，任時磊，高曉玲，桂志國

(1. 中北大學 生物醫學成像與影像大數據山西省重點實驗室，山西 太原 030051；2. 臨汾市人民醫院，山西 臨汾 041000)

0 引 言

數字圖像是生活中重要的信息來源，然而，圖像在產生、傳輸、使用的過程中，不可避免地會引入噪聲，由于噪聲往往與圖像高頻信息相關，因此，在保留重要特征的同時去除噪聲是一項值得研究的工作. 圖像去噪在圖像增強、圖像分割、目標識別等領域扮演著重要的角色[1-3]，如何在去除噪聲的同時保留圖像邊緣和紋理細節，一直是研究的重點. 在圖像空間域經典的去噪算法有雙邊濾波、偏微分方程[4]、非局部均值[5]、導向濾波[6]等. 其中，文獻[4]提出的各向異性擴散方法也稱為P-M模型，將像素值看作熱流，把圖像梯度引入擴散系數，在圖像邊緣處，擴散系數小，在平滑處，擴散系數大，從而有效保護了圖像邊緣，該模型自提出以來，就受到研究人員的廣泛關注.

近年來，國內外眾多研究學者基于P-M模型提出了很多改進算法. Ochotorena等[7]提出將各向異性擴散與加權導向濾波[8]結合，利用加權平均實現最大擴散，同時保留圖像的強邊緣，解決了細節光暈的問題. 方政等[9]提出基于多方向中值濾波的各向異性擴散算法，通過局部方差和圖像梯度改進擴散系數，結合中值濾波算法，在濾波過程中注重對邊緣細節的保持，提高了圖像質量. Chao等[10，11]提出在各向異性擴散模型中針對不同類型圖像，在模型中加入灰度方差，進一步保留了圖像的紋理區域. 董嬋嬋等[12]在P-M模型的基礎上，將差分曲率引入各向異性擴散方程中，較好區分了圖像邊緣、平坦區域和噪聲，自適應擴散取得了較好的降噪效果. Hossein等[13]分析圖像紋理的振蕩屬性，在殘差圖像中利用紋理信息構造紋理檢測算子，結合梯度模值改進擴散系數，有效保護了圖像紋理信息. 這些方法雖然有效改進了P-M模型，但處理結果仍存在受噪聲影響大、圖像細節易丟失等問題. 近年來，研究學者陸續提出新的去噪方法，非局部均值算法借鑒多幅相同圖像疊加取均值的方法可有效去除噪聲的原理，利用圖像中大量的冗余信息，在圖像搜索窗中尋找局部相似塊，根據圖像塊相似程度確定去噪權重，達到較好的去噪效果，但其較長的運行時間和邊緣易模糊的缺點會影響其大規模應用. 研究人員從各向異性擴散和非局部均值的優缺點出發，提出了很多改進算法. Yuan[14]提出在非局部均值模型中加入亮度、空間位置和梯度信息減弱P-M模型的階梯效應，有效保護了圖像細節. 陳強等[15]在擴散系數函數中，將像素的片相似代替圖像梯度，同時擴展到彩色圖像的去噪中，取得了較好效果，但依舊存在去噪不徹底、邊緣易模糊等問題. Yang等[16]提出的基于非局部均值理論的各向異性擴散(NL P-M)算法研究P-M模型的不穩定性，在迭代過程中，先將噪聲圖像使用非局部均值算法處理，然后將處理后圖像的梯度值代替噪聲圖像的梯度值，取得了較好的去噪效果，但其也存在過度平滑和效率低的問題. 為了提高非局部均值算法的效率，研究人員提出了很多解決方案，其中著名的有基于積分圖[17-19]的方法，其先計算出像素的積分圖，用加減運算代替圖像塊歐氏距離計算，以存儲空間換計算時間的方式，降低算法的復雜度.

本文改進算法受文獻[16]的啟發，在P-M模型中，將非局部均值引入迭代過程中，利用圖像中大量相似冗余信息有效地平滑圖像和去除噪聲，同時在擴散模型中，引入殘差圖像局部能量紋理檢測算子，更好地保護了圖像紋理和細節；針對非局部均值算法效率低的問題，采用積分圖的方法降低算法復雜度，減少了算法運行時間.

1 算法原理

1.1 P-M模型及其不穩定性分析

經典的非線性P-M模型根據圖像梯度模值自適應擴散，其迭代公式為

t∈(0,T),

(1)

(2)

g(|?u|)=e-(|?u|/k)2,

(3)

式中：k為梯度擴散閾值，k選擇較大值時，|?u|/k有較小值，g(|?u|)趨于1，圖像擴散趨于各向同性擴散，圖像模糊，相反，k選擇較小值時，g(|?u|)趨于0，抑制圖像擴散，但容易造成圖像去噪不徹底的問題.

文獻[16]研究P-M模型的不穩定性，將圖像分解為梯度方向η和切線方向ξ，得到

(4)

(5)

則有

(6)

(7)

式中：uη和uηη分別為圖像u在梯度方向的一階和二階偏導數；uxx，uyy，uxy分別為圖像u(x,y,t)的二階方向導數.類似得，u(x,y,t)在切線方向的二階偏導數為

(8)

將式(1)分解為

[g(|?u|)·ux]x+[g(|?u|)·uy]y=

g(|?u|)(uxx+uyy)=

g(|?u|)(uηη+uζξ)=

[g′(|?u|)|?u|+g(|?u|)]uηη+g(|?u|)uζξ,

(9)

式中：g(|?u|)永遠為正值，意味著正向擴散和圖像平滑.然而，g′(|?u|)|?u|+g(|?u|)可以是正值也可以是負值，導致不穩定的結果. 再加上噪聲的影響，各向異性擴散容易出現階梯效應等問題.

P-M模型擴散的離散形式為

un+1=un+Δtdiv[g(|?u|)?u]=

un+Δtg(|?u|)|?u|，

(10)

式中：Δt為時間步長，其范圍大小為0.05≤Δt≤0.25，?u梯度可用有限差分表示為

(11)

且有

g(|?u|)|?u|=g(|?Nui,j|)|?Nui,j|+

g(|?Sui,j|)|?Sui,j|+g(|?Wui,j|)|?Wui,j|+

g(|?Eui,j|)|?Eui,j|.

(12)

1.2 非局部均值算法及積分圖加速

非局部均值理論自提出以來，一直受到研究人員的廣泛關注，其根據多幅相同圖像疊加可有效去除圖像中隨機噪聲的原理，利用圖像中像素的鄰域相似性，確定圖像搜索窗內每個像素的去噪權值，計算公式為

(13)

式中：v為含噪聲圖像；uNL為去噪后的圖像；權值w(i,j)表示像素i和像素j之間的相似性，其由像素i和像素j為中心的鄰域像素決定，保證其相似性度量的準確性.

(14)

Z(x)為歸一化系數，計算公式為

(15)

式中：h為平滑參數，其值越大圖像越平滑，其中

(16)

i為圖像鄰域內索引.

假設圖像中共有N個像素，搜索窗大小為D×D，計算像素相似度的匹配窗口大小為d×d(d=2ds+1)，ds為匹配窗鄰域半徑，因此，計算像素鄰域相似度的時間復雜度為O(d2), 每個像素在搜索窗內共進行D2次計算，所以，非局部均值算法的時間復雜度為O(Nd2D2).由于非局部均值算法的時間復雜度較高，研究人員提出了積分圖加速的方法.

在圖像鄰域求相似度中，定義積分圖像為

(17)

其中：

(18)

兩個像素領域的相似性可以用下式計算，大大減少運行時間.

St(x1-ds-1,x2-ds-1)-St(x1+ds,x2-ds-1)-

St(x1-ds-1,x2+ds)).

(19)

因此，積分圖加速后，非局部均值算法的時間復雜度降為O(ND2).

2 本文算法

(20)

式(10)改進為

(21)

文獻[13]為自適應選擇擴散去噪算法，在各向異性擴散中，引入殘差局部能量紋理檢測算子，可以更準確地檢測圖像紋理和細節，同時降低階梯偽影. 將殘差圖像定義為非局部均值處理結果與再經過P-M模型擴散結果的差值圖像，如式(22)所示

(22)

PR=Ps+Pn,

(23)

式中：PR,Ps,Pn分別為uR,us,un的局部能量.因為Pn是近似恒定的，可以將其認為是PR的能量偏移強度，紋理信息Ps在PR中有較大的值.PR的計算公式為

PR(x,y)=

(24)

式中：μ(·)為均值算子；w(x,y)(x,y)為歸一化的高斯算子，經過高斯運算，殘差圖中的噪聲部分被平滑，紋理部分更加突出.

如圖 1 所示，在殘差局部能量圖中，紋理區域對應圖像中較大值，使紋理區域有了較好的提取，隨著迭代次數的增加，局部能量算子可以有效地濾除噪聲的影響，從而更準確地檢測紋理細節.

圖 1 局部能量及第一次和第二次迭代的殘差局部能量

由于梯度算子往往較大，為了使局部能量的取值更為合理，使局部能量歸一化到梯度范圍

(25)

式中：PRScaled為紋理檢測算子.

將紋理檢測算子與梯度算子結合，即為文獻[13]中提到的算法，在圖像擴散過程中，可以提升圖像保護紋理細節的能力，但容易出現去噪不徹底的問題. 本文結合上述兩種算法的優點，在NL P-M算法中，改進擴散函數，添加局部能量紋理檢測算子，提出改進的非局部均值各向異性擴散去噪算法，其基本流程為：

(26)

3) 根據式(24)計算殘差圖像的局部能量，并根據式(25)將局部能量歸一化到圖像梯度范圍，其結果如圖 1 所示.

(27)

5) 重復1)～4)，直到獲取最終的去噪結果.

3 實驗結果與分析

為驗證所提算法的有效性，使用大小為512×512的標準灰度圖像進行測試. 使用Matlab 2014a進行圖像實驗，計算機配置為：CPU為Intel(R) Core(TM) i7-4770 @3.50 GHz，內存為8 GB. 因為本文算法是P-M模型結合非局部均值算法的改進，所以采用經典的P-M模型算法、自適應選擇擴散去噪算法和NL P-M算法為對比算法.

3.1 實驗參數設置

本文算法參數較多，當面對不同類型圖像和不同強度噪聲時，需要適當手動調整以獲取最佳的去噪結果. 對于擴散閾值k，取較大值時，容易使去噪圖像平滑模糊，取較小值時，去噪不徹底. 當圖像噪聲強度較大時，應適當選擇較大的擴散閾值. 平滑參數h在非局部均值濾波中發揮作用，其值較大時，容易使圖像過度平滑，因此，當圖像含有較多紋理細節時，其值應適當減小. 紋理控制閾值β在文獻[13]中的范圍為β≥1，β值越大，保留的紋理信息越多.高斯核參數σg為文獻[13]中的參考值.在所提算法中，非局部均值中平滑參數h和紋理控制閾值β的值對算法有較大影響.

3.2 主觀視覺效果分析

本文在標準灰度測試圖像中加入均值為0，標準差為20的高斯白噪聲，實驗效果及局部放大圖像如圖 2 所示.

圖2 添加標準差σ=20時，不同圖像去噪結果

仔細觀察圖 2 及局部放大圖像可以看到，經典P-M算法在強邊緣處取得較好結果，在圖像平滑處，容易受到噪聲的影響，出現階梯效應和少量斑點偽影；文獻[13]所提算法中結合紋理檢測算子，弱邊緣和紋理處的效果都比較好，但在平滑處出現去噪不徹底的問題，視覺效果較差；NL P-M算法中，雖然視覺效果較好，但在紋理處出現過平滑的問題；本文算法有效克服了上述缺點，在平滑圖像和保留紋理細節方面都有較好的表現，在去除噪聲的同時保留了更多的紋理細節，比如在Lena和Barbara的放大圖像中更好保留的紋理，Airplane圖像中較為清晰的數字.

3.3 客觀評價指標分析

在圖像去噪領域，常用的客觀評價標準有峰值信噪比(PSNR)和平均結構相似度(MSSIM)[21].

PSNR的定義如下

(28)

式中：m和n分別為圖像行數和列數；u為恢復圖像；u0為原始噪聲圖像.PSNR值越大，代表恢復結果越好.

結合圖像亮度、對比度和結構信息，具有符合人眼視覺系統特性的結構相似度(SSIM)定義為

(29)

在實際評估中，常采用MSSIM來估計整幅圖像的結構相似度，即為

(30)

式中：μx和μy分別為x和y的均值；σx和σy為標準差；σxy為協方差；J為局部窗口的總個數.MSSIM值越大，表示去噪結果越接近原始圖像.

因主觀視覺評價受觀測者個體差異性較大，本節從PSNR和MSSIM兩個客觀指標來分析所提算法的去噪效果，結果如表1，表2 所示.

表1 添加高斯白噪聲，經不同算法處理的PSNR

從表1 和表2 可以看出，在標準測試圖像中添加標準差為20和30的高斯白噪聲，本文所提算法相對P-M算法、文獻[13]算法和NL P-M算法，都有最好的峰值信噪比和平均結構相似度，因此，從視覺效果和客觀評價指標都驗證了所提算法的有效性. 通過對不同噪聲強度的去噪處理和選擇最為常用的Lnea，Barbara，Airplane標準測試圖像，驗證了所提算法的有效性.

表2 添加高斯白噪聲后經不同算法處理的MSSIM

3.4 算法時間復雜度分析

基于各向異性擴散的方法具有運行速度快的優勢，但由于結合了非局部均值算法，使本文算法時間復雜度大大增加，因此，采用基于積分圖方法的算法加速. 設圖像共有N個像素，迭代次數為iter，下面分析對比算法與所提算法的時間復雜度.

P-M算法中用有限差分的方式求每個像素的梯度，并計算散度算子，因此，其時間復雜度為

T(n)=O(4iterN).

(31)

文獻[13]在擴散函數中加入紋理檢測算子，增加了紋理檢測算子的計算時間，其時間復雜度為

T(n)=O(4iterN+N).

(32)

NL P-M算法在每次迭代中都有非局部均值算法，根據前面的分析，非局部均值的時間復雜度為O(Nd2D2)，因此，NL P-M的時間復雜度為

T(n)=O(iterND2d2+4iterN).

(33)

本文算法的時間復雜度和NL P-M算法的時間復雜度基本相同，但是由于采用了積分圖加速方法，時間復雜度降低，此外，本算法和NL P-M算法迭代次數較少，算法耗時主要為非局部均值運行時長，時間復雜度為

T(n)=O(iterND2+4iterN).

(34)

所提算法結合矩陣操作，優化代碼后可進一步減少運行時間. 為了更直觀地觀察各種算法的運算效率，以大小為512×512的Lena圖像為例，將各種算法的運行時間記錄如表3 所示，可見本文所提算法相比NL P-M算法，運行時間大幅度減少.

表3 Lena圖像在標準差為20時，各種算法的運行時長

4 結 論

本文算法在P-M模型算法的基礎上，將非局部均值算法引入迭代過程中，即在各向異性擴散前進行預處理，減弱了P-M模型的不穩定性；同時在擴散系數函數中融合殘差圖像局部能量算子，在擴散過程中，更準確地保護了紋理細節，提出一種改進的非局部均值各向異性擴散去噪算法，并且在非局部算法中，運用積分圖的方法和使用矩陣運算代替傳統像素遍歷循環，極大減少了算法耗時，解決了NL P-M算法過度平滑和運行時間較長的問題. 從視覺效果和客觀評價指標可以看出，所提算法優于P-M模型算法、文獻[13]算法和NL P-M算法，在有效去除噪聲的同時，更多地保留了紋理細節信息. 但是，本文主要針對在圖像中最容易引入的高斯白噪聲，對于其他類型的噪聲去除，將是后繼研究的重點.