陳文斌，王燕玲，王 娜，劉 祎，張 權(quán)，邵燕靈

(1. 中北大學(xué) 生物醫(yī)學(xué)成像與影像大數(shù)據(jù)山西省重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051；2. 中北大學(xué)，理學(xué)院，山西 太原 030051)

0 引 言

數(shù)字圖像在采集、傳輸和儲(chǔ)存等過程中不可避免地引入了噪聲，而在眾多應(yīng)用領(lǐng)域都需要清晰的圖像，因此圖像復(fù)原具有重要的研究意義. 目前，常用的圖像去噪算法有： 小波濾波法，均值濾波法、高斯濾波法、中值濾波法、非局部均值濾波法、以及基于偏微分方程的圖像濾波法等. 在這些去噪算法中，基于偏微分方程(PDE)的圖像去噪算法是最具代表性的一類方法，如Koenderink[1]提出的熱擴(kuò)散方程和Perona和Malik[2]提出的各向異性擴(kuò)散PM模型等. 熱擴(kuò)散方程是一個(gè)各向同性擴(kuò)散方程，去噪的同時(shí)會(huì)引起圖像邊緣的模糊. 為了克服這個(gè)缺點(diǎn)，1990年P(guān)erona和Malik提出了各向異性擴(kuò)散PM模型，其利用圖像的局部特征(梯度模)來決定擴(kuò)散系數(shù)，去噪的同時(shí)很好地保持了邊緣. 但由于PM模型的解是不穩(wěn)定的，會(huì)產(chǎn)生階梯效應(yīng)，所以許多正則化技術(shù)被應(yīng)用于PM模型中. 如Catte等[3]提出了一個(gè)正則化PM模型，其解是連續(xù)依賴原始圖像的，它用高斯平滑后的圖像梯度模來決定擴(kuò)散系數(shù)，去噪效果明顯優(yōu)于原始PM模型. 但PM模型以及它的正則化模型都只受圖像梯度模的控制，然而圖像梯度模不能代表全部的圖像局部特征，所以Weickert[4]又引入了擴(kuò)散張量，求解加入擴(kuò)散張量的模型是一個(gè)適定問題. Gliboa等[5]提出了前向后向自適應(yīng)擴(kuò)散模型，可以在平滑噪聲的同時(shí)加強(qiáng)圖像的細(xì)節(jié)特征； 后來，他們又提出了保持斜坡的復(fù)擴(kuò)散模型[6]. 文獻(xiàn)[7]提出了一個(gè)加權(quán)曲率保持PDE模型，其權(quán)重由擴(kuò)散張量決定，該方法能有效地保持邊緣和曲率幾何結(jié)構(gòu). 文獻(xiàn)[8]中的擴(kuò)散系數(shù)含有殘余局部方差和梯度模，使得對(duì)圖像去噪時(shí)不僅有效地保持了邊緣結(jié)構(gòu)還保持了紋理和細(xì)節(jié)信息. PM模型的擴(kuò)散效果主要依賴于四個(gè)參數(shù)： 迭代次數(shù)、時(shí)間步長(zhǎng)、邊緣閾值參數(shù)以及擴(kuò)散系數(shù)函數(shù)，文獻(xiàn)[9-11]等做了很多優(yōu)化這些參數(shù)的工作.

二階PM模型在去噪時(shí)會(huì)產(chǎn)生階梯效應(yīng). 為了消除階梯效應(yīng)，很多學(xué)者提出了高階PDE去噪技術(shù)[12-16]. 例如You-Kaveh[13]提出了Y-K模型，該模型在邊緣切線方向和梯度方向進(jìn)行不同程度的擴(kuò)散，在梯度方向的擴(kuò)散系數(shù)比在邊緣切線方向的擴(kuò)散系數(shù)小. YK的改進(jìn)模型[17]中，用差分曲率和片相似模構(gòu)造了新的擴(kuò)散系數(shù)，去噪時(shí)不僅有效地減輕了斑點(diǎn)效應(yīng)和階梯效應(yīng)，還很好地保持了邊緣和細(xì)節(jié)信息. 另外一個(gè)經(jīng)典的四階模型是Lundervold和Tai[14]提出的LLT模型. 四階模型去除高頻噪聲的速率比二階的快，這意味著四階模型會(huì)過度平滑階躍性邊緣，除此之外，還會(huì)在結(jié)果圖像中產(chǎn)生斑點(diǎn)噪聲. 為了讓PM模型在去噪的同時(shí)能夠保持圖像的邊緣等細(xì)節(jié)信息，許多技術(shù)被運(yùn)用于PM模型中，如梯度矢量流(GVF)[18]、改進(jìn)的梯度矢量流(IGVF)[19].

鑒于在去噪的同時(shí)，原始PM模型雖然能較好地保持圖像邊緣，但不能保持細(xì)節(jié)信息，還產(chǎn)生階梯效應(yīng)等問題，本文構(gòu)造了一個(gè)結(jié)合局部方差[20]和梯度模的新擴(kuò)散系數(shù)，梯度模將邊緣和平坦區(qū)域以及細(xì)節(jié)區(qū)域區(qū)分，局部方差又將細(xì)節(jié)區(qū)域與平坦區(qū)域區(qū)分. 另外，本文所提模型用自適應(yīng)的邊緣閾值參數(shù)代替了PM模型中的常數(shù)邊緣閾值參數(shù)[21]，對(duì)圖像的不同區(qū)域進(jìn)行不同強(qiáng)度的擴(kuò)散. 在邊緣區(qū)域，差分曲率[22]大，邊緣閾值參數(shù)小，擴(kuò)散強(qiáng)度小，有利于保持邊緣； 在獨(dú)立噪聲點(diǎn)和平坦區(qū)域，差分曲率小，邊緣閾值參數(shù)大，擴(kuò)散強(qiáng)度大，有利于去噪.

1 各向異性擴(kuò)散PM模型及其改進(jìn)模型

1.1 熱擴(kuò)散方程

傳統(tǒng)的熱擴(kuò)散方程是圖像去噪的有效工具，去噪后的結(jié)果圖像相當(dāng)于原始噪聲圖像和高斯核做卷積，擴(kuò)散強(qiáng)度在圖像的每個(gè)方向是一樣的，沒有考慮圖像的局部特征. 熱擴(kuò)散方程為

(1)

式中：I0是原始噪聲圖像，div是散度算子，是梯度算子，c(常數(shù))是擴(kuò)散系數(shù).

熱擴(kuò)散方程有一些缺陷：

1)圖像的邊緣會(huì)模糊；

2)圖像的邊緣總是移動(dòng)，偏離“自然”位置.

1.2 各向異性擴(kuò)散PM模型及其改進(jìn)模型

由于熱擴(kuò)散方程對(duì)邊緣和非邊緣區(qū)域進(jìn)行同樣強(qiáng)度的平滑，不利于邊緣的保持. 在1990年，Perona和Malik提出了各向異性擴(kuò)散PM模型，PM算法為

(2)

式中： 擴(kuò)散系數(shù)c(|I|)是關(guān)于|I|的單調(diào)遞減非負(fù)函數(shù). Perona和Malik給定了兩種常用的擴(kuò)散系數(shù)為

c(|I|)

(3)

c(|I|)

(4)

式中：k是邊緣閾值參數(shù)，控制擴(kuò)散的程度.

PM模型中的擴(kuò)散系數(shù)c是自適應(yīng)的，在邊緣處，梯度模|I|大，c比較小，有利于保持邊緣； 在非邊緣區(qū)域處，梯度模|I|小，c比較大，有利于去除噪聲.

雖然PM模型在去噪的同時(shí)能保持邊緣信息，但是由于PM模型是病態(tài)的，擴(kuò)散過程不穩(wěn)定，會(huì)產(chǎn)生階梯效應(yīng).

為了消除PM模型的病態(tài)性質(zhì)，Catte提出了PM模型的正則化模型[3]，即在每一次迭代開始前先用一個(gè)高斯波器處理噪聲圖像. Catte模型為

(5)

式中：Gσ*I表示用一個(gè)標(biāo)準(zhǔn)差為σ的高斯核和圖像I做卷積. 這個(gè)方法可以有效地去除噪聲，但是選擇合適的σ比較困難.

后來，Chao和Tsai提出了一個(gè)可以保持細(xì)節(jié)的改進(jìn)各向異性擴(kuò)散模型，記為CTPM模型. 該模型的擴(kuò)散系數(shù)結(jié)合了局部方差和梯度模，可以在去噪的同時(shí)較好地保持邊緣和細(xì)節(jié)信息. CTPM模型為

(6)

(7)

2 本文改進(jìn)的各向異性擴(kuò)散模型

2.1 差分曲率

文獻(xiàn)[21]提出的差分曲率可以較好地區(qū)分邊緣、平坦區(qū)域以及獨(dú)立噪聲，其表達(dá)式為

D=‖Iηη|-|Iξξ‖,

(8)

式中：Iηη和Iξξ分別表示圖像在梯度方向和邊緣切線方向的二階導(dǎo)數(shù). 梯度和邊緣切線方向的方向向量分別為

(9)

從而

Iη=I·η=

(10)

(11)

同理可得

(12)

其中

Ix=I(x+1,y)-I(x,y),

Iy=I(x,y+1)-I(x,y),

Ixx=I(x+1,y)+I(x-1,y)-2·I(x,y),

Iyy=I(x,y+1)+I(x,y-1)-2·I(x,y),

在邊緣處，|Iηη|較大，|Iξξ|較小，所以差分曲率D較大； 在平坦區(qū)域處，|Iηη|和|Iξξ|都較小，所以差分曲率D較小； 在獨(dú)立噪聲點(diǎn)處，|Iηη|和|Iξξ|都較大，但是幾乎相等，所以差分曲率D較小. 因此根據(jù)差分曲率可以很好地在噪聲和平坦區(qū)域中將邊緣提取出來.

2.2 自適應(yīng)閾值函數(shù)k′

傳統(tǒng)的PM算法，在迭代過程中使用的邊緣閾值參數(shù)k是同一個(gè)常數(shù)(一般是直接給定，或者通過噪聲估計(jì)，如Canny算子). 我們知道參數(shù)k決定了結(jié)果圖像的平滑程度，當(dāng)k過小會(huì)導(dǎo)致降噪不充分，結(jié)果圖像中仍有未去除的噪聲； 當(dāng)k較大會(huì)造成圖像的過平滑，邊緣以及細(xì)節(jié)信息得不到保護(hù)，事實(shí)上根據(jù)經(jīng)驗(yàn)很難選取一個(gè)合適的k使算法對(duì)圖像降噪的同時(shí)充分保護(hù)邊緣以及細(xì)節(jié)等信息. 所以，在擴(kuò)散系數(shù)中設(shè)置一個(gè)自適應(yīng)的k是很有必要的. 在對(duì)圖像濾波過程中，我們的目的是在去噪的同時(shí)能保持邊緣和細(xì)節(jié)信息，也就是在平坦區(qū)域擴(kuò)散系數(shù)能自動(dòng)地增大，在邊緣和細(xì)節(jié)區(qū)域擴(kuò)散系數(shù)能自動(dòng)地減小.

差分曲率可以較好地判斷像素點(diǎn)是在平坦區(qū)域、邊緣區(qū)域還是孤立噪聲點(diǎn)，所以可以用差分曲率來自適應(yīng)地控制k，從而對(duì)不同區(qū)域的像素點(diǎn)采用不同的邊緣閾值參數(shù)k′，實(shí)行不同強(qiáng)度的擴(kuò)散，本文提出的自適應(yīng)邊緣閾值參數(shù)k′的表達(dá)式為

(13)

式中：k0是給定的常數(shù)，大于傳統(tǒng)PM模型中的k，D是歸一化差分曲率. 具體分析如下：

在第一次迭代開始之前，給定k0，在擴(kuò)散過程中根據(jù)圖像的局部特征，k′會(huì)自適應(yīng)地調(diào)整，在平坦區(qū)域和獨(dú)立噪聲點(diǎn)處，邊緣閾值參數(shù)k′較傳統(tǒng)PM模型的邊緣閾值參數(shù)k大，增大了擴(kuò)散強(qiáng)度，從而可以更好地去除噪聲； 在邊緣區(qū)域，邊緣閾值參數(shù)k′較傳統(tǒng)PM模型的邊緣閾值參數(shù)小，降低了擴(kuò)散強(qiáng)度，從而可以更好地保持邊緣等細(xì)節(jié).

2.3 新擴(kuò)散系數(shù)

傳統(tǒng)的PM模型利用梯度模來區(qū)分邊緣區(qū)域和平坦區(qū)域以及細(xì)節(jié)區(qū)域，但是梯度模卻不能有效地區(qū)分細(xì)節(jié)區(qū)域和平坦區(qū)域，去噪時(shí)丟失了一部分細(xì)節(jié)信息，本文中引用文獻(xiàn)[19]中的局部方差來有效地區(qū)分細(xì)節(jié)區(qū)域和平坦區(qū)域，使得在去噪的同時(shí)，不僅可以保持圖像的邊緣信息而且可以保持圖像的細(xì)節(jié)信息.

為了簡(jiǎn)單起見，把改進(jìn)的PM模型簡(jiǎn)記為IPM模型. IPM模型為

(14)

首先，定義一個(gè)大小為3×3的模板，則It(x,y)的局部方差為

(15)

歸一化局部方差為

(16)

下面通過實(shí)驗(yàn)驗(yàn)證IPM模型的效果.

3 實(shí)驗(yàn)結(jié)果與分析

為了驗(yàn)證本文提出的IPM算法的有效性和可行性，以大小為512×512添加高斯噪聲的Lena圖像和barbara圖像為例，進(jìn)行實(shí)驗(yàn)，并將IPM算法與PM算法和CTPM算法進(jìn)行對(duì)比分析.圖1(a) 為L(zhǎng)ena原圖像，圖1(b)為barbara原圖像.

圖 1 原圖圖像Fig.1 Original image

圖2(a)和圖2(b)為含噪標(biāo)準(zhǔn)差為10的Lena和barbara圖像，圖2(c)和圖2(d)為含噪標(biāo)準(zhǔn)差為15的Lena和barbara圖像. 本文及對(duì)比算法均在Windows 7旗艦版32位操作系統(tǒng)，Intel? Pentium? CPU G2020 @ 2.90 GHz，內(nèi)存為4.00 GB的電腦上進(jìn)行，編程使用的工具為MATLABR2012a.

圖 2 噪聲圖像Fig.2 Noise image

圖 3 為PM算法、CTPM算法和本文IPM算法處理含噪標(biāo)準(zhǔn)差為10的Lena圖像和barbara圖像的對(duì)比結(jié)果，其中圖3(a)～圖3(c)為三種算法對(duì)Lena圖像的處理結(jié)果圖，圖3(d)～圖3(f)為三種算法對(duì)barbara圖像的處理結(jié)果圖.

圖 3 噪聲標(biāo)準(zhǔn)差為10的各去噪算法對(duì)比圖Fig.3 Comparison of denoising algorithms for Lena image and barbara image (σ=10)

圖 4 對(duì)應(yīng)圖3中各子圖的放大圖Fig.4 Correspond to the enlargement of each subgraph in Fig.3

圖 4 表示圖3中各算法對(duì)應(yīng)的局部區(qū)域放大圖像，由圖4(a)可知，PM算法在降噪的時(shí)候會(huì)保持邊緣不被破壞，但帽檐位置的褶皺(細(xì)節(jié)信息)會(huì)變模糊； 由圖4(b) 可知，CTPM算法在擴(kuò)散系數(shù)中加入了梯度模和局部方差，實(shí)驗(yàn)證明CTPM算法性能優(yōu)于PM算法，保持邊緣的同時(shí)能夠有效地保持細(xì)節(jié)信息，帽檐位置處的褶皺比較明顯，但仍殘留部分噪聲. 從圖4(c)可知，IPM算法相比PM算法和CTPM算法不僅能夠有效地去除圖像中含有的噪聲，還可以很好地保持圖像邊緣以及帽檐等細(xì)節(jié)信息.

圖 5 為PM算法、CTPM算法和本文IPM算法處理含噪標(biāo)準(zhǔn)差為15的Lena圖像和barbara圖像的對(duì)比結(jié)果.

圖 5 噪聲標(biāo)準(zhǔn)差為15的各去噪算法對(duì)比圖Fig.5 Comparison of denoising algorithms for Lena image and barbara image at σ=15

其中圖5(a)～圖5(c)為三種算法對(duì)Lena圖像的處理結(jié)果圖，圖5(d)～圖5(f)為三種算法對(duì)barbara圖像的處理結(jié)果圖；圖6表示圖5中各算法對(duì)應(yīng)的局部區(qū)域放大圖像. 從圖6(d)可明顯看到，PM算法在去噪時(shí)會(huì)產(chǎn)生階梯效應(yīng)，形成虛假的邊緣，視覺效果不佳. 從圖6(e)中看到，CTPM算法相比PM算法而言，雖然適當(dāng)?shù)販p輕了階梯效應(yīng)，也保持了邊緣，但是去噪效果不徹底，仍有部分噪聲殘留. 從圖6(c)中看到，IPM算法不僅沒有產(chǎn)生明顯的階梯效應(yīng)，在很好地濾掉噪聲的同時(shí)，不僅保持了圖像的邊緣，還保持了細(xì)節(jié)信息. 因此，IPM算法相比PM算法和CTPM算法在噪聲去除、減輕階梯效應(yīng)和保持邊緣、細(xì)節(jié)信息方面等都有了較大的提升.

圖 6 對(duì)應(yīng)圖5中各子圖的放大圖Fig.6 Correspond to the enlargement of each subgraph in Fig.5

從上述分析可以看出，本文提出算法的去噪效果，即噪聲去除能力與邊緣、細(xì)節(jié)保持能力都明顯優(yōu)于其他兩種去噪算法，為了對(duì)算法進(jìn)行客觀、定量的評(píng)價(jià)，本文采用較常用的信噪比(Signal-to-Noise Ratio, SNR)、峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)和平均結(jié)構(gòu)相似度(Mean Structural Similarity Index，MSSIM)對(duì)各去噪算法進(jìn)行定量的描述. 這些參數(shù)的定義為

1) 信噪比(Signal-to-Noise Ratio, SNR)

式中：I是結(jié)果圖像；Id是原始無噪圖像.

SNR越高，結(jié)果圖像的質(zhì)量越好，噪聲越少.2) 峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)

PSNR在一定程度上反映了圖像處理前后數(shù)據(jù)變化情況的統(tǒng)計(jì)平均，其值越大，圖像失真越少.

3) 平均結(jié)構(gòu)相似度(Mean Structural Similarity Index，MSSIM)

SSIM(x,y)=

由表 1 分析可知，IPM去噪算法的信噪比(SNR)、峰值信噪比(PSNR)都高于其他對(duì)比算法，所以本文算法去噪效果較好，由于其平均結(jié)構(gòu)相似度(MSSIM)也高于其他對(duì)比算法，可知本文算法處理后的結(jié)果圖更接近真實(shí)圖像，有利于邊緣和細(xì)節(jié)等的保持. 通過主觀視覺觀察和質(zhì)量評(píng)價(jià)參數(shù)的定量評(píng)估，都驗(yàn)證了本文算法的可行性和優(yōu)越性.

表 1 各降噪算法的客觀評(píng)價(jià)Tab.1 Objective assessment of various denoising algorithms

4 結(jié)束語

傳統(tǒng)PM算法的擴(kuò)散系數(shù)僅依賴于圖像的梯度模，但梯度模不能很好地區(qū)分細(xì)節(jié)區(qū)域和平坦區(qū)域，在去噪的過程中會(huì)丟失細(xì)節(jié)信息，所以本文在傳統(tǒng)PM算法中加入了局部方差來區(qū)分細(xì)節(jié)區(qū)域和平坦區(qū)域，在去噪的過程中不僅能保持邊緣信息還能保持細(xì)節(jié)信息. 傳統(tǒng)PM模型中，整個(gè)擴(kuò)散過程的邊緣閾值參數(shù)是同一個(gè)固定常數(shù)，即對(duì)整個(gè)圖像進(jìn)行的是同等程度的擴(kuò)散，這顯然不太合理，人們希望，在平坦區(qū)域選擇比較大的，擴(kuò)散強(qiáng)度大一些，有利于去噪，在邊緣細(xì)節(jié)區(qū)域，擴(kuò)散強(qiáng)度比較小一些，有利于保持邊緣和細(xì)節(jié). 因此本文提出了一種自適應(yīng)邊緣閾值參數(shù)，通過差分曲率來自適應(yīng)地調(diào)節(jié)邊緣閾值參數(shù)，在邊緣處，差分曲率大，邊緣閾值參數(shù)小，擴(kuò)散強(qiáng)度小，有利于保持邊緣，在平坦區(qū)域和獨(dú)立噪聲點(diǎn)處，差分曲率小，邊緣閾值參數(shù)大，擴(kuò)散強(qiáng)度大，有利于去噪. 實(shí)驗(yàn)結(jié)果表明，本文的改進(jìn)IPM模型在去噪的同時(shí)，消除了階梯效應(yīng)，能很好地保持邊緣和細(xì)節(jié)等信息，無論在主觀視覺效果還是定量評(píng)價(jià)參數(shù)方面都是可行的.