王常虹，陳韜亦，屈楨深

(哈爾濱工業(yè)大學(xué) 空間控制與慣性技術(shù)研究中心，黑龍江 哈爾濱 150001)

醫(yī)學(xué)超聲圖像的去噪問(wèn)題一直是國(guó)內(nèi)外超聲成像技術(shù)的重要課題之一.由于其成像機(jī)制的限制，圖像清晰度不高是超聲成像的主要缺點(diǎn).此外，超聲圖像中特有的散斑噪聲不僅降低了超聲圖像的質(zhì)量，也使對(duì)圖像細(xì)節(jié)的識(shí)別與分析變得更加困難，因此，超聲圖像去噪一般要求有效抑制散斑噪聲，同時(shí)要保留對(duì)后期分析和診斷有用的圖像細(xì)節(jié)信息.超聲圖像去噪的主要難點(diǎn)在于:1)散斑噪聲可以大致看作為一種乘性噪聲［1-3］;2)噪聲的隨機(jī)性質(zhì)比較復(fù)雜;3)噪聲易與圖像細(xì)節(jié)相混，而圖像細(xì)節(jié)又復(fù)雜且多樣.

目前超聲圖像的去噪方法主要包括基于中值濾波、基于小波變換和基于擴(kuò)散方程等.基于中值濾波的方法根據(jù)圖像的局部統(tǒng)計(jì)特征來(lái)自動(dòng)選取濾波窗口內(nèi)各點(diǎn)權(quán)值及窗口的大小和形狀，盡管在保留圖像細(xì)節(jié)方面取得了一定效果，但對(duì)窗口的選擇很敏感，限制了處理效果.基于小波變換的方法將圖像變換到小波域，利用小波閾值處理將某些認(rèn)為是噪聲的系數(shù)丟棄，再逆變換回圖像，但如何選擇小波變換的尺度和閾值尚無(wú)確定方法.基于擴(kuò)散方程的去噪方法通過(guò)求解初始值為輸入圖像的非線性熱擴(kuò)散方程來(lái)實(shí)現(xiàn).在擴(kuò)散方程中，通過(guò)引入圖像特征，設(shè)計(jì)合適的擴(kuò)散系數(shù)來(lái)控制擴(kuò)散方程的擴(kuò)散行為，使得在平滑圖像的同時(shí)能夠保留甚至增強(qiáng)圖像的特征信息.

本文的基本思想是以超聲圖像為背景，借助多方向中值這種改進(jìn)中值濾波方式，結(jié)合由歸一化局部方差和圖像梯度組成的改進(jìn)擴(kuò)散系數(shù)，使濾波模型在去除散斑噪聲能力、保持邊緣能力以及迭代速度等指標(biāo)上比傳統(tǒng)各向異性擴(kuò)散算法有更好的效果.

1 現(xiàn)有算法分析與評(píng)價(jià)

1.1 傳統(tǒng)去噪算法

傳統(tǒng)的圖像平滑算法如均值濾波、中值濾波和高斯濾波等，沒(méi)有考慮圖像像素的形狀特征，即無(wú)論像素處于特征區(qū)域還非特征區(qū)域，都做同樣的濾波處理，其平滑結(jié)果等價(jià)于傳導(dǎo)系數(shù)為常量的熱擴(kuò)散方程，屬于各向同性擴(kuò)散，因而在去除噪聲的同時(shí)會(huì)模糊甚至破壞圖像的邊緣信息.

針對(duì)超聲圖像和合成孔徑雷達(dá)(synthetic aperture radar，SAR)圖像中特有的散斑噪聲，現(xiàn)在廣泛采用的濾波器有Lee濾波器和Frost濾波器等.這2種算法均為根據(jù)局部均值和方差的局部統(tǒng)計(jì)濾波算法.當(dāng)局部的數(shù)據(jù)比較均勻時(shí)，因?yàn)榇藭r(shí)的數(shù)據(jù)只含噪聲和變化較緩的信號(hào)，進(jìn)行有力的濾波.另一方面，當(dāng)局部的數(shù)據(jù)中存在較大變化時(shí)，僅進(jìn)行較輕微濾波或者不濾波，因?yàn)榇藭r(shí)數(shù)據(jù)中包含邊緣或者其他結(jié)構(gòu)變化.這種濾波策略對(duì)含噪聲的邊緣處理得不夠.

1.2 PM模型

與熱擴(kuò)散模型相比較，各向異性擴(kuò)散模型實(shí)際上是一個(gè)非線性的偏微分方程，由圖像的梯度來(lái)決定擴(kuò)散速度，能夠同時(shí)兼顧噪聲消除和特征保持2個(gè)方面.目前，以Perona-Malik(PM)模型為代表的各向異性平滑方法己在邊緣檢測(cè)、圖像增強(qiáng)、圖像分割以及目標(biāo)識(shí)別等領(lǐng)域得到了廣泛的應(yīng)用.

PM模型提出的各向異性擴(kuò)散方程為［4］

式中，x和y代表圖像I0的空間坐標(biāo)，div是散度算子，▽為梯度算子，t代表運(yùn)算時(shí)間，c(x)為擴(kuò)散系數(shù).PM模型提出了2類擴(kuò)散系數(shù):

式中，k稱為梯度閾值，一般根據(jù)實(shí)際圖像估計(jì)為常數(shù).在這里，邊緣采用常用的梯度微分算子來(lái)識(shí)別，可以理解為PM模型將邊緣檢測(cè)和噪聲去除很好地統(tǒng)一到了變分方法的偏微分方程中.

PM模型擴(kuò)散系數(shù)的選取應(yīng)符合如下基本原則:

1)平滑量的控制:在圖像特征多的區(qū)域平滑量應(yīng)該盡可能地小甚至不平滑;在圖像特征少的區(qū)域或沒(méi)有圖像特征的區(qū)域平滑量要大.

2)平滑方向的控制:沿圖像特征方向的擴(kuò)散量大，穿越圖像特征的擴(kuò)散量小甚至不平滑.

所以，各向異性擴(kuò)散系數(shù)c(x)是各向異性擴(kuò)散算法的關(guān)鍵所在，假設(shè)c(x)是一個(gè)常數(shù)，那么擴(kuò)散就是各向同性的，它的效果等同于線性低通濾波.

PM模型存在不足之處，集中表現(xiàn)為以下幾點(diǎn):

1)對(duì)于孤立噪聲點(diǎn)，PM模型的平滑效果較差.因?yàn)檫@類噪聲存在較大的梯度，將作為邊緣保留下來(lái).

2)擴(kuò)散系數(shù)的設(shè)計(jì)使PM模型的邊緣保持效果不理想.由于c(x)永遠(yuǎn)無(wú)法取值為0，即使是微小擴(kuò)散經(jīng)多次迭代后也會(huì)被放大，最終造成邊緣的模糊乃至損壞，而且所提供的2類擴(kuò)散系數(shù)都只是啟發(fā)性的，還沒(méi)有足夠的理論來(lái)支持其正確性和合理性.

3)PM模型簡(jiǎn)單地采用圖像梯度的單調(diào)遞減函數(shù)作為擴(kuò)散系數(shù)，擴(kuò)散系數(shù)中的梯度閾值k決定了擴(kuò)散系數(shù)的特性，很多情況下梯度閾值被看成常數(shù).然而，由于噪聲圖像中梯度具有很大的不確定性，并且隨著圖像平滑程度的增加，相應(yīng)的梯度不斷下降，所以k遞減才能有效地保持邊緣，因?yàn)檫吘壉３值臈l件是|▽I|＞k.

4)從數(shù)學(xué)角度來(lái)看，各向異性擴(kuò)散方程本身在數(shù)學(xué)上是病態(tài)的，不能保證解的存在唯一性.

1.3 Catté模型

為解決PM模型的不足，Catté等［5］指出在噪聲較強(qiáng)的圖像中，梯度的準(zhǔn)確計(jì)算將變得非常困難，在一些強(qiáng)噪聲區(qū)域也會(huì)有較大的梯度，此時(shí)噪聲會(huì)被當(dāng)成圖像的邊界被保持下來(lái)，從而產(chǎn)生偽邊緣，通過(guò)高斯平滑方式來(lái)計(jì)算梯度信息能更好地估計(jì)局部梯度，提高抗噪性，并證明了這種改進(jìn)為方程解的存在、正則化和唯一性提供了充分條件.然而超聲圖像的散斑噪聲是空間相關(guān)且非對(duì)稱分布的，空間不變且各向同性的高斯平滑過(guò)程背離了各向異性擴(kuò)散算法的本質(zhì)，而且高斯平滑過(guò)程還會(huì)導(dǎo)致圖像結(jié)構(gòu)偏離原始位置.再者，高斯卷積核標(biāo)準(zhǔn)差σ和梯度閾值k的大小如何確定，Catté等沒(méi)有給出解決方案.

1.4 SRAD模型

Yu等［6］在討論了各向異性擴(kuò)散方程和Lee濾波、Frost濾波之間聯(lián)系的基礎(chǔ)上，將異性擴(kuò)散方程應(yīng)用到散斑抑制處理中，提出了去除散斑噪聲的各向異性擴(kuò)散模型(speckle reducing anisotropic diffusion，SRAD).Yu等修正了擴(kuò)散系數(shù)，使擴(kuò)散方程能夠根據(jù)圖像噪聲的情況而調(diào)整擴(kuò)散系數(shù)，并且能夠?qū)D像的細(xì)節(jié)信息更加敏感.

但SRAD模型也存在明顯的缺陷:模型中尺度函數(shù)q0(t)是由初始圖像中盡可能大的均勻區(qū)域計(jì)算得到，模型的關(guān)鍵是如何選取圖像中一個(gè)盡可能大的、合適的均勻區(qū)域.這一區(qū)域的選取往往很大程度上影響著擴(kuò)散結(jié)果，給實(shí)驗(yàn)結(jié)果帶有較大的偶然性.此外，SRAD模型所使用的局部統(tǒng)計(jì)信息實(shí)際上是各向同性的，這也背離了各向異性擴(kuò)散算法的本質(zhì).

1.5 關(guān)于參數(shù)k的最佳估計(jì)

Voci等［7］給出了2種估計(jì)參數(shù)k的方法:基于形態(tài)學(xué)的方法和基于P－范數(shù)的方法，從尋求參數(shù)k的最佳估計(jì)方法這點(diǎn)入手，來(lái)改進(jìn)各向異性擴(kuò)散系數(shù).基于形態(tài)學(xué)的k估計(jì)方法由于在每次迭代過(guò)程中要對(duì)整幅圖像做開、閉操作，計(jì)算成本較大.因?yàn)镻－范數(shù)有單調(diào)下降這一特性，所以P－范數(shù)估計(jì)能保證k遞減，且計(jì)算量小，缺點(diǎn)是平滑時(shí)模糊程度大，細(xì)節(jié)保持能力較差.以上方法存在問(wèn)題的原因是它們僅僅采用梯度來(lái)度量空間細(xì)節(jié)，得到的局部細(xì)節(jié)描述并不完整，因而細(xì)節(jié)保持能力不理想.

2 提出算法原理與實(shí)現(xiàn)

綜合以上各種去噪算法的優(yōu)缺點(diǎn)，著重分析以下2點(diǎn):

1)使用哪種方式進(jìn)行去噪前的預(yù)處理，使散斑噪聲對(duì)后續(xù)梯度檢測(cè)及局部統(tǒng)計(jì)信息提取的影響盡可能小;

2)使用哪種自適應(yīng)手段，能夠克服PM算法中梯度閾值為常數(shù)以及SRAD算法中需要預(yù)先手工指定均勻區(qū)域等缺陷帶來(lái)的魯棒性差的問(wèn)題.文獻(xiàn)［8］中，仍然使用傳統(tǒng)PM模型中k為常數(shù)的Tukey biweight等［9］擴(kuò)散系數(shù)方程和標(biāo)準(zhǔn)中值濾波作為一次迭代的全過(guò)程，對(duì)處理低信噪比的分子圖像顯示了較好的平滑效果，但是k的魯棒性差等問(wèn)題依然沒(méi)有解決.受其啟發(fā)，考慮到將中值濾波引入到各向異性擴(kuò)散中的良好效果和應(yīng)用前景，確定使用對(duì)超聲圖像散斑抑制效果好的改進(jìn)中值濾波和對(duì)局部細(xì)節(jié)敏感的改進(jìn)各向異性擴(kuò)散方程相結(jié)合的算法方案.

2.1 多方向中值濾波

根據(jù)Catté的說(shuō)明，對(duì)噪聲污染的圖像進(jìn)行PM擴(kuò)散前應(yīng)先進(jìn)行高斯濾波，然而高斯濾波破壞了各向異性擴(kuò)散的本質(zhì)，且會(huì)導(dǎo)致圖像結(jié)構(gòu)偏離原始位置.而同屬各向同性擴(kuò)散濾波的均值濾波在平滑圖像的同時(shí)，對(duì)細(xì)節(jié)和邊緣的模糊作用也非常嚴(yán)重，以上方法均為線性濾波器.于是選擇保留細(xì)節(jié)能力較為突出的非線性濾波器中的典型代表，即中值濾波，因?yàn)閷?duì)于非對(duì)稱分布的噪聲，中值濾波的性能要優(yōu)于線性濾波器，而且中值濾波保持了圖像結(jié)構(gòu)的位置.

考查中值濾波的各種改進(jìn)形式，發(fā)現(xiàn)方向中值濾波［10］對(duì)邊緣方向特別敏感，保留細(xì)節(jié)能力強(qiáng).該算法提取所有一維中值濾波器中的最大輸出作為濾波結(jié)果，減少了邊緣模糊的程度，但是其濾波結(jié)果易受噪聲，特別是高密度噪聲的影響.

為了減少噪聲對(duì)中值濾波器的影響，同時(shí)考慮到邊緣不同于背景，其灰度分布存在明顯的方向性，綜合以上思路提出多方向中值濾波算法.該方法的主要思想是針對(duì)不同方向的中值濾波器，采用不同尺寸的局部窗口，結(jié)合各方向的濾波結(jié)果來(lái)估計(jì)最終的濾波器輸出.

為簡(jiǎn)單起見(jiàn)，只考慮水平和垂直邊緣2種情況.假設(shè)水平中值濾波器的窗口大小是Mx×Nx(Mx＜Nx)，垂直中值濾波器的窗口大小是My×Ny(My＜Ny)，故水平方向的中值濾波器的作用范圍是一個(gè)扁的長(zhǎng)方形，垂直方向的中值濾波器的作用范圍是一個(gè)細(xì)高的長(zhǎng)方形，這樣保證了該濾波器盡可能多地包含了同樣灰度性質(zhì)的邊緣點(diǎn)，減少了中值濾波器的模糊作用.例如點(diǎn)P是水平邊緣上的一點(diǎn)，其水平中值濾波器的輸出應(yīng)該比垂直中值濾波器的輸出更接近真實(shí)值，因?yàn)樗椒较蛏吓cP點(diǎn)相同的點(diǎn)要比垂直方向的多.如果選擇水平中值濾波器的輸出作為該點(diǎn)的濾波輸出，則去噪的同時(shí)邊緣得以保持.對(duì)背景點(diǎn)，其水平中值濾波器的輸出與垂直中值濾波器的輸出的差別不大，濾波結(jié)果為2個(gè)中值濾波器輸出的平均值.在實(shí)際使用中，由于不知道邊緣的方向，假設(shè)P點(diǎn)本身雖可能受噪聲影響，但仍比較接近真值，這樣就可以選擇與P點(diǎn)灰度值接近的濾波器的輸出作為濾波結(jié)果.當(dāng)然也使本方法有了一定的局限，如果P點(diǎn)受噪聲的影響過(guò)大，而其所在的邊緣又很細(xì)，該點(diǎn)的濾波效果不太好，可以通過(guò)考慮P點(diǎn)鄰域的邊緣方向加以克服.

設(shè)含噪圖像為I(x，y)，濾波結(jié)果為g(x，y)，算法流程如下:

1)計(jì)算水平方向和垂直方向的中值濾波結(jié)果h(x，y)和v(x，y);

2)分別計(jì)算點(diǎn)(x，y)的灰度值I(x，y)與h(x，y)和v(x，y)的差:kx=|I(x，y)－h(huán)(x，y)|，ky= |I(x，y)－v(x，y)|，然后分情況考慮:

①如果kx＜ky，若水平鄰點(diǎn)(x，y－1)屬于水平方向的邊緣，那么g(x，y)=h(x，y)，此時(shí)點(diǎn)(x，y)屬于水平方向的邊緣，否則由2個(gè)濾波器結(jié)果的線性插值更新g(x，y)=ah(x，y)+bv(x，y)，其中a=ky/(kx+ky)，b=kx/(kx+ky)，點(diǎn)(x，y)可以是水平或垂直邊緣.如果點(diǎn)(x，y)位于圖像邊上，沒(méi)有水平鄰點(diǎn)，則直接由h(x，y)更新.

②如果kx＞ky，使用1)中相同的方式對(duì)點(diǎn)(x，y)的垂直鄰點(diǎn)(x－1，y)進(jìn)行考查.

③若kx=ky，那么g(x，y)=0.5h(x，y)+0.5· v(x，y)，點(diǎn)(x，y)可以是水平或垂直邊緣.

多方向中值濾波算法在去除噪聲的同時(shí)，由于考慮了邊緣的方向性，提升了中值濾波算法的邊緣保持能力.雖然該算法只計(jì)算水平和垂直2個(gè)方向，在實(shí)際應(yīng)用中完全可以以相同方式考慮更多方向的濾波結(jié)果，從而提高算法的性能.相對(duì)于中值濾波算法，該方法相當(dāng)于同時(shí)進(jìn)行多個(gè)中值濾波，其算法的復(fù)雜度并沒(méi)有大的增加，算法結(jié)構(gòu)也完全可以應(yīng)用原中值濾波的結(jié)構(gòu)，可以滿足超聲圖像的實(shí)時(shí)性要求.

2.2 歸一化局部均值和梯度

設(shè)I為原圖像，以任意一點(diǎn)I(x，y)為中心，以(2m+1，2n+1)為大小的矩形區(qū)域內(nèi)的局部灰度均值為

局部方差為

圖像沒(méi)有噪聲污染時(shí)，在平滑區(qū)域內(nèi)像素點(diǎn)的局部均值σ2(x，y)接近像素灰度值，從而局部方差很小;而在邊緣區(qū)域內(nèi)，像素灰度值變化活躍，部分像素點(diǎn)灰度值與局部均值之間的差值較大，使局部方差σ2(x，y)很大.但如果有噪聲存在，平滑區(qū)域和邊緣區(qū)域的像素點(diǎn)局部方差σ2(x，y)都變大，此時(shí)局部方差就不能有效地度量空間細(xì)節(jié)，于是采用歸一化局部方差(x，y)來(lái)代替 σ2(x，y)，并用(x，y)|▽I(x，y)|來(lái)度量空間細(xì)節(jié).其中:

對(duì)擴(kuò)散系數(shù)進(jìn)行如下改進(jìn):

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

首先給出本算法的離散化數(shù)值迭代形式，考慮四鄰域情況［4］:

其中為了保證迭代穩(wěn)定性，0≤λ≤0.25.

擴(kuò)散系數(shù)和歸一化局部方差在每次迭代時(shí)均更新一次，其中:

以及

為驗(yàn)證算法效果，在Matlab7.0上對(duì)本文算法進(jìn)行了仿真實(shí)驗(yàn).實(shí)驗(yàn)計(jì)算機(jī)采用AMD Athlon 64× 2處理器，主頻1.70 GHz，內(nèi)存1 G.分別使用人工圖像，現(xiàn)實(shí)圖像和臨床超聲圖像進(jìn)行對(duì)比驗(yàn)證.Lee濾波器和Frost濾波器本質(zhì)上均為各向同性擴(kuò)散，擴(kuò)散時(shí)對(duì)邊緣的保持能力差;Catté模型和Voci模型由于其固有的缺陷，在濾除散斑噪聲時(shí)效果并不令人滿意.為突出對(duì)比性，對(duì)PM模型分別取k=20和k=40，同時(shí)選取去除散斑噪聲效果良好的SRAD模型一起和本文提出算法模型進(jìn)行對(duì)比.

為取得相對(duì)最優(yōu)效果和迭代穩(wěn)定性，參照文獻(xiàn)［4，6］的參數(shù)設(shè)置，SRAD中時(shí)間步長(zhǎng)尺度Δt= 0.05，PM中時(shí)間步長(zhǎng)尺度Δt=1/7.本文算法中多方向中值濾波取垂直和水平2個(gè)方向，窗口大小分別為5×1和1×5，局部方差計(jì)算時(shí)使用3×3窗口，各向異性擴(kuò)散中λ=0.25.

3.1 人工圖像

如圖1所示，對(duì)一幅大小為320×320的人工圖像應(yīng)用4個(gè)去噪模型.首先對(duì)圖1(a)加入更具有一般意義的零均值，方差為0.3的乘性高斯白噪聲，迭代次數(shù)均為10次.圖1(c)～(f)分別顯示了4個(gè)模型的濾波結(jié)果.通過(guò)與其他3個(gè)模型對(duì)比，本文算法模型在去除乘性噪聲和邊緣保持方面得到了更好的效果.可以看出k值的選取能夠直觀地反映在迭代速度上，圖1(d)比圖1(c)有更好的去噪效果，但是對(duì)于產(chǎn)生的大方差點(diǎn)的濾除效果并不好.SRAD的邊緣保持能力得到了初步驗(yàn)證，可以從各種形狀的邊緣保持情況看出來(lái).但是由于迭代次數(shù)的限制，SRAD對(duì)乘性噪聲濾除的優(yōu)勢(shì)并沒(méi)有顯現(xiàn)出來(lái).相比而言，在同樣的迭代次數(shù)下，本文算法在視覺(jué)上得到了最好的去噪效果，雖然有些局部邊緣受到了過(guò)平滑的影響而顯得模糊.

圖1 人工圖像去噪效果對(duì)比Fig.1 Comparison of noise denoising capability for synthetic image

為了定量評(píng)估算法性能，從峰值信噪比(peak signal noise ratio，PSNR)和均方誤差(mean square error，MSE)2個(gè)方面對(duì)4個(gè)模型進(jìn)行比較.使用Multi_var表示加入乘性白噪聲的方差，且方差范圍在0.2～0.4［10］，結(jié)果如表1、2所示.PSNR和MSE的定義如下:

式中:X(i，j)是原始無(wú)噪圖像的像素值，Y(i，j)是加噪以后圖像的像素值，Z(i，j)是濾波后圖像的像素值.

表1 乘性高斯噪聲下的PSNRTable 1 PSNR Result for multiplicative Gaussian noise corrupted images

表2 乘性高斯噪聲下的MSETable 2 MSE Result for multiplicative Gaussian noise corrupted images

由表1、2可得，在迭代次數(shù)相同的條件下，從PSNR和MSE來(lái)看，本文算法都得到了最好的結(jié)果，說(shuō)明與原始圖像相似度最高.k=40時(shí)比k=20時(shí)效果要好，而SRAD由于迭代次數(shù)的限制，優(yōu)勢(shì)并沒(méi)有顯現(xiàn)，效果最差.

3.2 現(xiàn)實(shí)圖像

Loupas［11］指出對(duì)數(shù)壓縮超聲圖像(即超聲顯示設(shè)備圖像)斑點(diǎn)噪聲是與信號(hào)相關(guān)的噪聲，提出顯式對(duì)數(shù)壓縮超聲圖像仿真模型:

式中，f0為原信號(hào)，f為觀測(cè)信號(hào)，n為零均值標(biāo)準(zhǔn)方差為σn的高斯白噪聲.本節(jié)基于式(11)作為散斑噪聲模型.此外，為了比較不同方法間的邊緣保持能力，使用FOM(figure of merit)定義［12］如下:

式中，Nreal和Nideal分別是原始無(wú)噪圖像和加噪濾波圖像經(jīng)過(guò)邊緣探測(cè)之后的邊緣點(diǎn)數(shù)，di是第i個(gè)濾波圖像邊緣點(diǎn)和與之最接近的原始圖像邊緣點(diǎn)之間的歐氏距離，e是常數(shù)，取值為1/9.FOM的變化范圍是［0，1］，當(dāng)2幅圖完全一致時(shí)取最大值.

使用經(jīng)典的大小為256×256的Lena灰度圖像來(lái)進(jìn)行進(jìn)一步的性能比較.對(duì)圖2(a)加入方差為0.05的散斑噪聲，如圖2(b)所示.圖2(c)～(f)為4個(gè)模型的濾波結(jié)果，迭代次數(shù)分別為20、10、200、3次.圖2(g)～(k)為使用Canny邊緣探測(cè)器對(duì)圖2(b)～(f)取得的邊緣，參照文獻(xiàn)［6］，Canny邊緣探測(cè)器中高斯核標(biāo)準(zhǔn)差為0.1.

圖2 Lena圖像去噪效果對(duì)比Fig.2 Comparison of noise denoising capability for standard image Lena

由于迭代次數(shù)都能使模型達(dá)到最優(yōu)的FOM值，故圖2(c)～(f)和圖2(h)～(k)分別反映了各個(gè)模型真實(shí)的去噪效果和邊緣保持能力.

如表3所示，在大小均為512×512的4幅標(biāo)準(zhǔn)灰度圖像進(jìn)行比較，給出了各個(gè)模型濾波結(jié)果的PSNR最大值和對(duì)應(yīng)的迭代次數(shù).從PSNR最大值看，本文算法得到的濾波結(jié)果與原始圖像最為接近，另外3個(gè)模型的最佳濾波結(jié)果對(duì)應(yīng)的PSNR略小于本文算法，濾波質(zhì)量稍差.從迭代次數(shù)來(lái)看，本文算法次數(shù)最少，SRAD模型的次數(shù)最多.

表3 加噪圖像PSNR和迭代次數(shù)的比較Table 3 Comparison of PSNR and iteration times for speckled images

比較這4個(gè)模型迭代一次所需要的時(shí)間.以Mandrill(512×512)為例，PM(k=20)、PM(k=40)、SRAD和本文算法所需要的時(shí)間分別為0.295 7、0.294 4、0.187 9和1.083 s.所以綜合單次迭代時(shí)間和最佳濾波效果對(duì)應(yīng)的迭代次數(shù)，PM(k=40)模型能在最短的時(shí)間達(dá)到最優(yōu)濾波效果，本文算法次之，SRAD最慢.

以256×256的Lena圖像為例，圖3顯示了在受不同方差條件的散斑噪聲污染情況下，4個(gè)模型所能達(dá)到的最佳邊緣保持效果所對(duì)應(yīng)的FOM值.可以看到本文算法濾波后的邊緣效果要優(yōu)于其他3個(gè)模型.

圖3 不同噪聲條件下FOM的比較Fig.3 FOM comparison under variable noise condition

3.3 臨床超聲圖像

圖4、5分別是實(shí)際采集的普通腫瘤和乳腺腫瘤超聲圖像，實(shí)際超聲圖像中主要的噪聲形式是散斑噪聲，散斑噪聲是一種乘性噪聲，并且局部空間相關(guān)且非對(duì)稱分布.PM(k=20)、PM(k=40)、SRAD和本文算法的迭代次數(shù)分別為20、15、300和10次.

圖4 普通腫瘤超聲圖像去噪效果對(duì)比Fig.4 Comparison of noise denoising capability for common tumor ultrasound image

圖5 乳腺腫瘤超聲圖像去噪效果對(duì)比Fig.5 Comparison of noise denoising capability for breast tumor ultrasound image

從圖4、5中可以看出4個(gè)模型都對(duì)散斑噪聲有抑制作用，效果更好的是本文算法.PM模型算法保留了比較多的噪聲殘留，SRAD濾波后圖像細(xì)節(jié)呈“階梯”狀分布，腫瘤組織邊緣保留能力有限.

4 結(jié)論

本文得出結(jié)論如下:

1)多方向中值濾波算法在去除噪聲的同時(shí)，由于考慮了邊緣的方向性，提升了中值濾波算法的邊緣保持能力;

2)在改進(jìn)各向異性擴(kuò)散模型中，使用歸一化局部方差和圖像梯度組成的擴(kuò)散系數(shù)有效地提高了空間細(xì)節(jié)描述的準(zhǔn)確程度;

3)多組仿真實(shí)驗(yàn)表明本文提出的基于多方向中值濾波和改進(jìn)各向異性擴(kuò)散的去噪算法比傳統(tǒng)PM和SRAD模型有更好的濾除超聲圖像噪聲和邊緣保持能力.

此外，如何在擴(kuò)散過(guò)程中根據(jù)平滑后圖像的質(zhì)量，設(shè)計(jì)一個(gè)好的算法迭代停止準(zhǔn)則，是今后進(jìn)一步的研究方向.

［1］SALINAS H M，F(xiàn)ERNANDEZ D C.Comparison of PDE-based nonlinear diffusion approaches for image enhancement and denoising in optical coherence tomography［J］.IEEE Trans on Medical Imaging，2007，26(6):761-771.

［2］TAY，P C，GARSON C D，ACTON S T，et al.Ultrasound despeckling for contrast enhancement［J］.IEEE Trans on Image Processing，2010，19(7):1847-1860.

［3］RABBANI H，VAFADUST M，ABOLMAESUMI P，et al.Speckle noise reduction of medical ultrasound images in complex wavelet domain using mixture priors［J］.IEEE Trans on Biomedical Engineering，2008，55(9):2152-2160.

［4］PERONA P，MALIK J.Scale-space and edge detection using anisotropic diffusion［J］.IEEE Trans on Pattern Analysis and Machine Intelligence，1990，12(7):629-639.

［5］CATTE T，LIONS P，MOREL J，et al.Image selective smoothing and edge detection by nonlinear diffusion［J］.SIAM J Numerical Anal，1992，29:182-193.

［6］YU Y，ACTON S.Speckle reducing anisotropic diffusion［J］.IEEE Trans on Image Processing，2002，11(11): 1260-1270.

［7］VOCI F，EIHO S，SUGIMOTO N，et al.Estimating the gradient threshold in the Perona-Malik equation［J］.IEEE Signal Processing Magazine，2004，21(3):39-46，65.

［8］LING J，BOVIK A C.Smoothing low-SNR molecular images via anisotropic median-diffusion［J］.IEEE Trans on Medical Image，2002，21(4):377-384.

［9］BLACK M J，SAPIRO G，MARIMONT D H，et al.Robust anisotropic diffusion［J］.IEEE Trans on Image Processing，1998，7(3):421-432.

［10］CZERWINSKI R N，JONES D L，WILLIAM D O.Ultrasound speckle reduction by directional median filtering［C］//IEEE International Conference on Image Processing.Washington DC，USA，1995:358-361.

［11］LOUPAS T，MCHICKEN W N，ALLAN P L.An adaptive weighted median filter for speckle suppression in medical ultrasonic image［J］.IEEE Trans on Circuits＆Systems，1989，36(1):129-135.

［12］RANJANI J J，THIRUVENGADAM S J.Dual-tree complex wavelet transform based SAR despeckling using interscale dependence［J］.IEEE Trans on Geoscience and Remote Sensing，2010，48(6):2723-2731.