陳文山

（福建衛(wèi)生職業(yè)技術(shù)學(xué)院，福州 350101）

0 引言

因?yàn)槌暢上窬哂邪踩?、?shí)時性和廉價性的特點(diǎn)，在臨床診斷中廣泛應(yīng)用。在超聲診斷中希望通過三維分割[1]、三維重建[2]，從而更全面獲取目標(biāo)腔體整體信息。由于超聲圖像存在獨(dú)有的斑點(diǎn)噪聲，為提高三維分割效率，希望在盡可能保護(hù)邊緣的前提下對超聲圖像進(jìn)行濾波。

目前經(jīng)典的圖像增強(qiáng)濾波方法基本上可分為兩大類：空域處理法和頻域處理法。然而，線性濾波能夠有效去除噪聲，但同時會使圖像的邊界模糊；中值濾波能有效去除脈沖噪聲，但因?yàn)榇翱趦?nèi)的排序處理破壞了空間相鄰信息和結(jié)構(gòu)，所以在去除噪聲的同時使圖像細(xì)節(jié)丟失、棱角模糊?；诟黜?xiàng)異性擴(kuò)散方程的濾波方法可以去除噪聲、邊界保留，但由于各向異性擴(kuò)散實(shí)質(zhì)上為迭代過程，要通過多次迭代，計算量大，運(yùn)算速度慢。

考慮到后續(xù)在超聲圖像三維重建的應(yīng)用需要，希望尋求一種有效去除噪聲的同時又不丟失邊緣信息，并且具有較高運(yùn)算速度的超聲圖像濾波算法。本文在各向異性擴(kuò)散方程的基礎(chǔ)上，結(jié)合中值濾波高效濾除斑點(diǎn)噪聲，各向同性擴(kuò)散方程在圖像同質(zhì)區(qū)域高速擴(kuò)散的特點(diǎn)，根據(jù)圖像梯度，設(shè)計一種基于邊緣保持的醫(yī)學(xué)超聲圖像濾波方法。

1 算法原理及實(shí)現(xiàn)方法

1.1 擴(kuò)散方程模型

圖像處理中的各向異性擴(kuò)散理論源于對流體擴(kuò)散分析的類比研究[3]。它將圖像灰度值I看作流體中的溶質(zhì)，根據(jù)熱擴(kuò)散理論，為了使?jié)舛雀叩娜苜|(zhì)和其周圍區(qū)域達(dá)到一種平衡，溶質(zhì)會向周圍區(qū)域擴(kuò)散，擴(kuò)散方程為在這里，D為擴(kuò)散張量（即擴(kuò)散系數(shù)），?t為對時間的導(dǎo)數(shù)。各向異性擴(kuò)散方程實(shí)質(zhì)為具有各項(xiàng)異性擴(kuò)散系數(shù)的偏微分方程，而對各向異性擴(kuò)散方程的迭代過程就是求解初始值為輸入圖像的偏微分方程。

1.1.1 各向同性擴(kuò)散模型（常系數(shù)熱傳導(dǎo)方程模型）

Marr等人將原圖像u0（x，y）與幾個帶不同參數(shù)σ的高斯核Gσ（x，y）作卷積，得到該圖像在若干尺度下的平滑圖像u（x，y）。各向同性擴(kuò)散模型在圖像的平滑過程中具有兩個特性：1）當(dāng)尺度由小變大時，圖像的平滑效果增強(qiáng)，圖像的細(xì)節(jié)特征逐漸消失；2）在大尺度上保留的圖像特征，其空間位置會展寬并發(fā)生漂移，使圖像的邊界、端點(diǎn)模糊。

1.1.2 各向異性擴(kuò)散模型

由于各向同性擴(kuò)散模型在平滑圖像的同時也破壞了圖像的一些重要的信息特征。為了克服其缺點(diǎn)，Perona和Malik給出了各向異性擴(kuò)散模型。與各向同性擴(kuò)散方程模型相比較，它實(shí)際上是一非線性拋物型微分方程。因其擴(kuò)散速度隨著梯度模的不同而改變，所以該模型為各向異性擴(kuò)散。然而，當(dāng)圖像含有的噪聲梯度較大時，異性擴(kuò)散模型對噪聲無法去除，帶有噪聲的邊界也會被破壞。

1.2 基于邊緣保持的醫(yī)學(xué)超聲圖像濾波方法

針對醫(yī)學(xué)超聲圖像噪聲大且在不同區(qū)域內(nèi)有不同的噪聲模型，而用于三維分割的圖像數(shù)據(jù)量大、對邊緣保持要求高的特點(diǎn)[4]，本文以圖像像素點(diǎn)梯度Y為閾值，采用綜合濾波方法。當(dāng)像素點(diǎn)梯度Y較大時，對圖像進(jìn)行中值去噪，去除高梯度噪聲；當(dāng)像素點(diǎn)梯度Y較小時，對圖像進(jìn)行各向異性擴(kuò)散，去除噪聲同時保持邊界；其余部分，對圖像進(jìn)行各向同性擴(kuò)散，去除隨機(jī)噪聲，提高濾波速度；

綜合濾波方法具體模型如下：

其中，Δu是各向同性的拉普拉斯算子；?u是?u(x,y,t))的簡寫，代表梯度算子；div是散度算子；是起擴(kuò)散濾波作用的退化擴(kuò)散項(xiàng)；u是當(dāng)前迭代的圖像，I是原始圖像，(1-g)(u-I)是保持初始圖像不過度擴(kuò)散的平衡控制項(xiàng)；g是用來檢測邊界和控制擴(kuò)散速度的尺度函數(shù)。Gσ是高斯函數(shù)，在此作為卷積算子對圖像進(jìn)行平滑濾波。

1.2.1 g函數(shù)的去噪原理

g函數(shù)在模型中起著增強(qiáng)邊界并控制擴(kuò)散速度的作用。當(dāng)圖像灰度差小時，梯度?(Gσ*u(x,y,t))表現(xiàn)得小，此時g函數(shù)接近于1，退化擴(kuò)散項(xiàng)充分發(fā)揮作用，擴(kuò)散幅度加大，噪聲被平滑。當(dāng)圖像灰度差大時，梯度?(GGσ*u(x,y,t))表現(xiàn)得大，此時g函數(shù)接近于0，退化擴(kuò)散項(xiàng)被控制，邊緣得到增強(qiáng)。

1.2.2 退化擴(kuò)散項(xiàng)的去噪原理

退化擴(kuò)散項(xiàng)是本文模型的核心，它有選擇的在圖像邊緣的切線上對圖像進(jìn)行擴(kuò)散濾波，在消除圖像噪聲的同時又增強(qiáng)了圖像的邊界。擴(kuò)散濾波是在邊緣的切線方向和梯度方向上的異性擴(kuò)散，因而對邊緣的保護(hù)作用有限；而退化擴(kuò)散濾波只在邊緣的切線方向上擴(kuò)散，而在梯度方向上不擴(kuò)散，從而對邊緣起到保護(hù)乃至增強(qiáng)的作用。

1.2.3 邊緣保持因子k的取值

由g函數(shù)的定義，知g函數(shù)由圖像高斯卷積后的梯度和常數(shù)k共同決定。前者由當(dāng)前處理的圖像決定，后者需要人為設(shè)定?？紤]到g函數(shù)的作用主要是限制退化擴(kuò)散項(xiàng)的擴(kuò)散行為，防止高梯度的圖像邊界被擴(kuò)散項(xiàng)破壞，而對于較低梯度的圖像同質(zhì)區(qū)，擴(kuò)散項(xiàng)的濾波功能應(yīng)該充分發(fā)揮。實(shí)驗(yàn)證明，對于帶有大量噪聲的醫(yī)學(xué)超聲圖像，g函數(shù)能夠更徹底的濾除噪聲，同時減少迭代次數(shù)從而節(jié)約了圖像處理時間。

1.2.4 平衡控制項(xiàng)原理

控制項(xiàng)具有使擴(kuò)散結(jié)果保持與原圖像重要信息一致的功能。它在合適的程度上降低了擴(kuò)散濾波對圖像信息的負(fù)面影響，但效果不甚理想。控制項(xiàng)和擴(kuò)散項(xiàng)之間的平衡由式子(1-g)控制。隨著(1-g)值的變化，模型會在圖像的不同位置呈現(xiàn)不同的效果：在同質(zhì)區(qū)域內(nèi)部，較小，g～1，模型有較強(qiáng)的擴(kuò)散作用；在尖銳處，例如邊緣，g～0，(1-g)～1，平滑后的圖像u會和原始圖像I保持一致。λ值在此處有調(diào)整平衡控制項(xiàng)功能的作用：λ值越大，原始圖像信息保持得越完整，但濾波效果越有損耗?？紤]到效率問題，本論文中一般取λ=1。

1.2.5 中值濾波及閾值選擇

本文模型中每次退化擴(kuò)散濾波前對圖像進(jìn)行的高斯濾波具有一定的平滑噪聲的作用。然而實(shí)踐中，僅使用“退化擴(kuò)散項(xiàng)＋平衡控制項(xiàng)”存在兩個問題：1）高斯函數(shù)參數(shù)σ不易控制，模型不穩(wěn)定；2）模型對邊緣選擇性增強(qiáng)，迭代完成后圖像邊界出現(xiàn)因過度增強(qiáng)產(chǎn)生的亮斑，在圖像同質(zhì)區(qū)域會產(chǎn)生新的細(xì)小散斑。本文設(shè)想引入中值濾波濾除隨機(jī)散粒噪聲。散粒噪聲點(diǎn)處局部梯度模具有較大值。當(dāng)局部梯度的模大于某個閾值時，在該點(diǎn)采用中值濾波；對于大部分局部梯度模小于閾值的像素點(diǎn)，依然采用各向異性擴(kuò)散。這樣不但去除噪聲，還能保持邊緣。實(shí)驗(yàn)中，圖像處理前由于大量噪聲的影響，代表圖像平滑的同質(zhì)區(qū)域小，代表噪聲的中等梯度區(qū)域大；圖像處理后噪聲被平滑，代表噪聲的中等梯度概率變小，代表圖像平滑的同質(zhì)區(qū)域大量增加，代表邊界的高梯度區(qū)依然保留。當(dāng)閾值過大時，中值濾波的作用無法發(fā)揮；當(dāng)閾值過小時，對圖像邊界有模糊作用且會有新的噪聲產(chǎn)生。經(jīng)過反復(fù)實(shí)驗(yàn)，本文取經(jīng)驗(yàn)值Ymax=120作為閾值。實(shí)驗(yàn)結(jié)果表明，引入中值濾波后，由于中值濾波對模型降噪作用的增強(qiáng)，處理圖像的迭代次數(shù)有所下降，縮短了圖像處理時間，提高了處理速度。

1.2.6 各向同性擴(kuò)散及閾值選擇

擴(kuò)散方程擴(kuò)散濾波中提到了拉普拉斯算子因?yàn)樗鼜膱D像邊緣的切線方向和梯度方向同時等大地進(jìn)行擴(kuò)散，所以被稱為各向同性擴(kuò)散。各向同性擴(kuò)散不能保留邊界，但對同質(zhì)區(qū)域內(nèi)的低梯度隨機(jī)噪聲有良好效果，更關(guān)鍵的是因?yàn)槠溥\(yùn)算相對簡單，它可以有效提高圖像處理的速度。隨著迭代次數(shù)的增加，低梯度的圖像區(qū)域明顯增加，而圖像梯度的變換也主要集中在低梯度區(qū)。這便從理論上說明了各向同性擴(kuò)散的閾值選擇也可以固定為一個很小的常數(shù)。實(shí)驗(yàn)證明，Ymin在0.5到1之間取值時，不但可以提高模型的速度，還可以配合中值濾波有效地改善因?yàn)橥嘶瘮U(kuò)散造成地邊緣亮斑和隨機(jī)散斑。當(dāng)Ymin取值過小時，各向同性濾波地效果不明顯；當(dāng)Ymin取值過大時，圖像會產(chǎn)生霧狀噪聲。

1.2.7 迭代終止準(zhǔn)則

方程迭代終止時間T的確定對擴(kuò)散方程濾波的結(jié)果有很大影響：如果T過小，則迭代的結(jié)果對初值圖像有很大的依賴性，使得圖像的很多噪聲不能被濾除；如果T過大，則迭代的結(jié)果就會被擴(kuò)散方程中固有的動力學(xué)特征所主宰，不但浪費(fèi)了時間，還會對圖像產(chǎn)生負(fù)效果。所以一個高效低計算量的迭代終止準(zhǔn)則對于模型的計算機(jī)實(shí)現(xiàn)是相當(dāng)必要的。當(dāng)圖像迭代不完全時，前后圖像間的變化大，當(dāng)圖像迭代基本完成后，前后圖像間的變化會極小。本文設(shè)想根據(jù)前后圖像間變化定義迭代終止準(zhǔn)則ISC，當(dāng)ISC滿足小于一定閾值時，圖像自適應(yīng)地停止迭代。

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

為了對本算法的有效性進(jìn)行驗(yàn)證，實(shí)驗(yàn)中選取了大量對比度較低，視覺效果較差[5]的醫(yī)學(xué)超聲圖像。圖1為其中一幅左心室聲像圖。

圖1 心臟超聲圖像實(shí)驗(yàn)結(jié)果

圖1(a)是原圖像，圖1(b)是增強(qiáng)Lee濾波器處理的圖像，圖1(c)是自適應(yīng)權(quán)值調(diào)節(jié)濾波器處理的圖像，圖1(d)是擴(kuò)散濾波處理的圖像，圖1(e)是各項(xiàng)異性擴(kuò)散濾波處理的圖像，圖1(f)是本文模型處理的圖像。通過對比可知，圖1(b)、圖1(c)、圖1(d)的濾波方法[6]，雖然一定程度上濾除了噪聲，但是整個圖像的清晰度有了明顯下降，輪廓也變得比較模糊。圖1(e)各項(xiàng)異性擴(kuò)散濾波處理后的圖像效果良好，既濾除了噪聲邊緣也得到了很好的保留。圖1(f)本文模型在達(dá)到上述效果的基礎(chǔ)上，還對邊緣起到了一定的銳化作用，本文模型在超聲圖像的去噪、保護(hù)邊緣上具有獨(dú)特的優(yōu)勢，不但能夠快速有效地消除噪聲，還能有效增強(qiáng)邊緣，同時邊緣的保持也與原始圖像最為吻合，邊緣線條更加平滑自然。邊緣信息的保護(hù)對超聲醫(yī)學(xué)圖像診斷具有十分重要的意義。

3 結(jié)論

大量實(shí)驗(yàn)結(jié)果表明，考慮圖像的梯度分布特征后，改進(jìn)的基于邊緣保持的醫(yī)學(xué)超聲圖像濾波方法，較現(xiàn)有的圖像濾波方法在濾除超聲圖像時具有更好地保護(hù)圖像邊緣、保留圖像細(xì)節(jié)信息的特點(diǎn)，是一種有效的醫(yī)學(xué)超聲圖像濾波方法。

[1]陳文山.基于曲面形變模型的心臟超聲圖像三維分割[J].福建電腦.2011(9):96-97.

[2]陳文山.基于三角網(wǎng)格模型的心臟超聲圖像三維表面重建[J].福建電腦.2014(11):20-22.

[3]Jemila Rose R,Allwin S. Denoising of Ultrasound Cervix Image Using Improved Anisotropic Diffusion Filter.[J].The West Indian medical journal,2015,64(4).

[4]陳文山.基于局部分形維數(shù)特征的醫(yī)學(xué)超聲圖像輪廓提取[J].福建電腦.2014(10):24-25.

[5]陳文山,汪天富,林江莉.基于相似度測量的醫(yī)學(xué)超聲圖像對比度增強(qiáng)[J].中國醫(yī)學(xué)影像技術(shù).2006(9):1432-1434.

[6]Qingzheng Wang,Shuai Li,Hong Qin,Aimin Hao. Robust multimodal medical image fusion via anisotropic heat diffusion guided low-rank structural analysis[J].Information Fusion,2015,26.