使用中值－各向異性擴散的超聲圖像去噪算法

王常虹，陳韜亦，屈楨深

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

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

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

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

1 現(xiàn)有算法分析與評價

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

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

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

1.2 PM模型

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

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

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

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

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

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

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

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

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

1)對于孤立噪聲點，PM模型的平滑效果較差.因為這類噪聲存在較大的梯度，將作為邊緣保留下來.

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

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

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

1.3 Catté模型

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

1.4 SRAD模型

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

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

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

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

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

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

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

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

2.1 多方向中值濾波

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

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

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

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

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

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

2)分別計算點(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，若水平鄰點(x，y－1)屬于水平方向的邊緣，那么g(x，y)=h(x，y)，此時點(x，y)屬于水平方向的邊緣，否則由2個濾波器結(jié)果的線性插值更新g(x，y)=ah(x，y)+bv(x，y)，其中a=ky/(kx+ky)，b=kx/(kx+ky)，點(x，y)可以是水平或垂直邊緣.如果點(x，y)位于圖像邊上，沒有水平鄰點，則直接由h(x，y)更新.

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

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

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

2.2 歸一化局部均值和梯度

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

局部方差為

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

對擴散系數(shù)進行如下改進:

3 仿真實驗及結(jié)果分析

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

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

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

以及

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

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

3.1 人工圖像

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

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

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

式中:X(i，j)是原始無噪圖像的像素值，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來看，本文算法都得到了最好的結(jié)果，說明與原始圖像相似度最高.k=40時比k=20時效果要好，而SRAD由于迭代次數(shù)的限制，優(yōu)勢并沒有顯現(xiàn)，效果最差.

3.2 現(xiàn)實圖像

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

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

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

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

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

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

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

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

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

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

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

3.3 臨床超聲圖像

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

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

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

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

4 結(jié)論

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

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

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

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

此外，如何在擴散過程中根據(jù)平滑后圖像的質(zhì)量，設(shè)計一個好的算法迭代停止準(zhǔ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.