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

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

0 引 言

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

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

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

1 算法原理

1.1 P-M模型及其不穩(wěn)定性分析

經(jīng)典的非線性P-M模型根據(jù)圖像梯度模值自適應(yīng)擴(kuò)散，其迭代公式為

t∈(0,T),

(1)

(2)

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

(3)

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

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

(4)

(5)

則有

(6)

(7)

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

(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|)永遠(yuǎn)為正值，意味著正向擴(kuò)散和圖像平滑.然而，g′(|?u|)|?u|+g(|?u|)可以是正值也可以是負(fù)值，導(dǎo)致不穩(wěn)定的結(jié)果. 再加上噪聲的影響，各向異性擴(kuò)散容易出現(xiàn)階梯效應(yīng)等問題.

P-M模型擴(kuò)散的離散形式為

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 非局部均值算法及積分圖加速

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

(13)

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

(14)

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

(15)

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

(16)

i為圖像鄰域內(nèi)索引.

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

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

(17)

其中：

(18)

兩個像素領(lǐng)域的相似性可以用下式計算，大大減少運(yùn)行時間.

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

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

(19)

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

2 本文算法

(20)

式(10)改進(jìn)為

(21)

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

(22)

PR=Ps+Pn,

(23)

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

PR(x,y)=

(24)

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

如圖 1 所示，在殘差局部能量圖中，紋理區(qū)域?qū)?yīng)圖像中較大值，使紋理區(qū)域有了較好的提取，隨著迭代次數(shù)的增加，局部能量算子可以有效地濾除噪聲的影響，從而更準(zhǔn)確地檢測紋理細(xì)節(jié).

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

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

(25)

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

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

(26)

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

(27)

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

3 實驗結(jié)果與分析

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

3.1 實驗參數(shù)設(shè)置

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

3.2 主觀視覺效果分析

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

圖2 添加標(biāo)準(zhǔn)差σ=20時，不同圖像去噪結(jié)果

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

3.3 客觀評價指標(biāo)分析

在圖像去噪領(lǐng)域，常用的客觀評價標(biāo)準(zhǔn)有峰值信噪比(PSNR)和平均結(jié)構(gòu)相似度(MSSIM)[21].

PSNR的定義如下

(28)

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

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

(29)

在實際評估中，常采用MSSIM來估計整幅圖像的結(jié)構(gòu)相似度，即為

(30)

式中：μx和μy分別為x和y的均值；σx和σy為標(biāo)準(zhǔn)差；σxy為協(xié)方差；J為局部窗口的總個數(shù).MSSIM值越大，表示去噪結(jié)果越接近原始圖像.

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

表1 添加高斯白噪聲，經(jīng)不同算法處理的PSNR

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

表2 添加高斯白噪聲后經(jīng)不同算法處理的MSSIM

3.4 算法時間復(fù)雜度分析

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

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

T(n)=O(4iterN).

(31)

文獻(xiàn)[13]在擴(kuò)散函數(shù)中加入紋理檢測算子，增加了紋理檢測算子的計算時間，其時間復(fù)雜度為

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

(32)

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

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

(33)

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

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

(34)

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

表3 Lena圖像在標(biāo)準(zhǔn)差為20時，各種算法的運(yùn)行時長

4 結(jié) 論

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