田玲 周盛 王曉春 葉青盛 王延群

300192天津，中國(guó)醫(yī)學(xué)科學(xué)院北京協(xié)和醫(yī)學(xué)院生物醫(yī)學(xué)工程研究所

眼科高頻超聲成像中斑點(diǎn)噪聲抑制算法的研究

田玲 周盛 王曉春 葉青盛 王延群

300192天津，中國(guó)醫(yī)學(xué)科學(xué)院北京協(xié)和醫(yī)學(xué)院生物醫(yī)學(xué)工程研究所

目的斑點(diǎn)噪聲是超聲圖像中存在的固有問題，而在眼科高頻超聲這種更為精細(xì)的超聲檢查中，有效地抑制斑點(diǎn)噪聲能提高圖像的質(zhì)量，有助于臨床醫(yī)生對(duì)病情的判別。方法提出了一種新的基于拉普拉斯（Laplacian）金字塔的多尺度斑點(diǎn)去噪方法。采用Laplacian金字塔，從斑點(diǎn)噪聲中分離出臨床圖像特征，根據(jù)每層子帶圖像不同尺度及特點(diǎn)，從小尺度到大尺度，首先采用改進(jìn)后的八方向各向異性斑點(diǎn)去噪（SRAD）去除圖像斑點(diǎn)，然后增強(qiáng)圖像的邊緣、細(xì)節(jié)及對(duì)比度等方面。該方法與傳統(tǒng)的SRAD濾波及相干增強(qiáng)濾波（CEDIF）進(jìn)行對(duì)比，采用等效視數(shù)及算法的時(shí)間耗費(fèi)對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行量化評(píng)估。結(jié)果與傳統(tǒng)SRAD濾波及CEDIF濾波方法相比，基于Laplacian金字塔的多尺度各向異性斑點(diǎn)去噪方法均高于前兩種方法（1.172 3 vs 1.122 3、0．929 3及0.864 0 vs 1.396 0、1．468 3）。結(jié)論本研究提出的基于Laplacian金字塔的多尺度各向異性斑點(diǎn)去噪方法在更有效地去除圖像斑點(diǎn)噪聲的同時(shí)，能很好地保存圖像邊緣及圖像細(xì)節(jié)等。

斑點(diǎn)噪聲；Laplacian金字塔；多尺度分析；各向異性濾波

Fund program：National Science&Technology Pillar Program(2012BAI13B02);the Central Level,Scientific Research Institutes for Basic R&D Special Fund Business

0 引言

超聲生物顯微鏡（ultrasound biomicroscope，UBM）是眼科超聲診斷儀中頻率最高的，也是目前超聲診斷設(shè)備中超聲頻率最高的設(shè)備之一，一般頻率在50 MHz以上。然而，本質(zhì)上UBM依然是采用了B超的成像原理[1-2]。UBM對(duì)眼前節(jié)的結(jié)構(gòu)和病變能很好地成像，可清晰地觀察角膜、虹膜、前后房角、睫狀體及懸韌帶等組織的結(jié)構(gòu)，進(jìn)而能診斷角膜病變等眼部疾病。然而由于超聲圖像固有斑點(diǎn)噪聲的影響，會(huì)削弱圖像對(duì)比度和細(xì)節(jié)分辨率，從而影響臨床判斷。因此，對(duì)含噪的超聲圖像進(jìn)行進(jìn)一步處理以得到更加清晰的圖像對(duì)于臨床診斷治療十分重要。

超聲圖像斑點(diǎn)去噪包括傳統(tǒng)的復(fù)合成像方法及常用的空域?yàn)V波[3]和頻域?yàn)V波[4]等。復(fù)合成像法[5]占用了系統(tǒng)較多的計(jì)算資源，目前通常只在高檔機(jī)上提供該功能。其他常用空域?yàn)V波分為單尺度、多尺度濾波等。然而這些方法雖然抑制了斑點(diǎn)噪聲，但卻損傷了包含圖像細(xì)節(jié)信息，如邊緣及醫(yī)療圖像中的病灶區(qū)等。因此，最大化地去除影響圖像的各種加性及乘性噪聲的同時(shí)，盡可能地保存圖像本身的紋理、邊緣特征及細(xì)節(jié)信息等重要特征信息不被影響很重要。

筆者提出了一種基于Laplacian金字塔[6]及改進(jìn)后的各向異性斑點(diǎn)去噪（speckle reduction anisotropic diffusion，SRAD）[7]的聯(lián)合去噪算法，從多個(gè)尺度對(duì)圖像進(jìn)行去噪、邊緣增強(qiáng)等。該方法的主要處理過(guò)程是對(duì)圖像進(jìn)行拉普拉斯金字塔分解，然后根據(jù)每層不同的特點(diǎn)，相應(yīng)地分別進(jìn)行斑點(diǎn)噪聲抑制或邊緣增強(qiáng)等，最后再將處理過(guò)的圖層進(jìn)行重建。

1 方法

1.1 Laplacian金字塔模型

Laplacian金字塔主要由分解和重構(gòu)兩部分組成，其中關(guān)鍵步驟分別記為Reduce和Expand。Reduce操作是對(duì)圖像進(jìn)行高斯濾波后，進(jìn)行一個(gè)在x和y方向上均為2的降采樣，每一層圖像縮小為上一層圖像的1/4，得到的這一組圖像為高斯金字塔；將得到的金字塔再進(jìn)行x和y方向?yàn)?的插值計(jì)算，濾波后與插值后的圖像相減，得到的即為殘差金字塔。殘差金字塔即包含了圖像在精細(xì)尺度及由插值得來(lái)的大尺度圖像之間的信息；而Expand則是進(jìn)行與之相反的過(guò)程。

假設(shè)輸入圖像為I，則高斯金字塔和Laplacian金字塔分別為G和L。

式中：G0為高斯金字塔第1層（原始圖像），G（l）為高斯金字塔的第l層，L（l）為L(zhǎng)aplacian金字塔的第l層。圖像的重構(gòu)即顛倒上述所敘步驟即可。

圖1為原始圖像進(jìn)行Laplacian分解之后的圖像（其原始圖像為UBM圖像）。

圖像分解為L(zhǎng)aplacian金字塔圖層后，根據(jù)斑點(diǎn)噪聲及圖像中某些元素的不同的頻率特性，其主要分布的圖層也不相同。斑點(diǎn)的噪聲頻率較高，通過(guò)觀察圖像證實(shí)它主要存在于低圖層里，即斑點(diǎn)噪聲主要存在于第1層；而第2層圖像，主要包含一些結(jié)構(gòu)元素；第3層，斑點(diǎn)噪聲基本上不存在了，此層中斑點(diǎn)噪聲的影響基本可忽略，可觀察到一些結(jié)構(gòu)信息（圖1）。根據(jù)每層殘差圖像的不同特點(diǎn)，需分別采取不同的濾波方法。

圖1 拉普拉斯金字塔

1.2 斑點(diǎn)噪聲抑制

由圖1可知，絕大多數(shù)斑點(diǎn)噪聲均被留在第1層圖像中。為了有效抑制斑點(diǎn)噪聲，同時(shí)保留邊緣的細(xì)節(jié)信息，筆者采用改進(jìn)后的八方向SRAD。傳統(tǒng)的SRAD[8]模型是基于PM[9]模型改進(jìn)的，其擴(kuò)散算法公式為

式中：I0（x，y）為原始圖像，▽表示梯度，div表示擴(kuò)散操作符，c（q）為擴(kuò)散系數(shù)（為非負(fù)單調(diào)遞減函數(shù)），?Ω代表Ω的邊界，n為?Ω邊線標(biāo)準(zhǔn)。其離散形式表現(xiàn)為

或

式中：q（x，y；t）為瞬時(shí)自變量系數(shù)，其表現(xiàn)形式為

q0（x，y）為尺度函數(shù)，其表達(dá)式為

式中：▽t為時(shí)間步，ρ為指數(shù)衰減率，q′0預(yù)先給定。

對(duì)于離散化的圖像，計(jì)算SRAD偏微分為

可得到圖像的更新函數(shù)為

在以上公式中，擴(kuò)散系數(shù)為SRAD模型的關(guān)鍵，其表示圖像處于圖像平坦的區(qū)域（c（q）→1）或細(xì)節(jié)邊緣區(qū)域（c（q）→0）。當(dāng)圖像處于均勻平坦的區(qū)域時(shí)，q→q0，SRAD表現(xiàn)為各向同性平滑濾波器；當(dāng)處于細(xì)節(jié)邊緣區(qū)域時(shí)，q?q0，SRAD表現(xiàn)為保持細(xì)節(jié)的功能。

筆者采用的改進(jìn)后的八方向SRAD，是在傳統(tǒng)SRAD擴(kuò)散方向上進(jìn)行了優(yōu)化，增加了左上、左下，右上、右下4個(gè)方向。原來(lái)上、下、左、右4個(gè)方向賦予的權(quán)重為2，左上、左下，右上、右下4個(gè)方向賦予的權(quán)重為1。由于擴(kuò)展方向以及權(quán)重的變化，使瞬時(shí)自變量系數(shù)中的數(shù)值發(fā)生改變，如下所示

在對(duì)第1層進(jìn)行斑點(diǎn)噪聲抑制后，為補(bǔ)償?shù)?層圖像中組織的線性結(jié)構(gòu)的缺損，采用了基于結(jié)構(gòu)張量[10]的相干增強(qiáng)濾波（coherence-enhancing diffusion, CEDIF）。還采用了對(duì)稱半正定擴(kuò)散張量來(lái)代替?zhèn)鹘y(tǒng)的標(biāo)量擴(kuò)散函數(shù)。擴(kuò)散張量D∈R2×2如下

式中：Ix、Iy分別為在x、y方向上的擴(kuò)散率。張量結(jié)構(gòu)的相角決定了擴(kuò)散的方向，而在該方向上的擴(kuò)散率是由局部相干信號(hào)決定的。

而在圖像的第3層中，斑點(diǎn)噪聲幾乎已經(jīng)完全可以忽略不計(jì)，可由圖1D（第3層殘差圖像）看出，它顯示圖像的整體結(jié)構(gòu)，但是分辨率相對(duì)比較低。因此在這一層采用非線性同態(tài)濾波器[11]，主要目的在于增強(qiáng)低對(duì)比度及增強(qiáng)微弱信號(hào)。同態(tài)濾波函數(shù)為

式中：常量c為步進(jìn)控制參數(shù)，β為高頻分量，α為低頻分量（濾波器由β和α控制改變），ε（ω）為傅里葉頻譜圖上角點(diǎn)到中間點(diǎn)的距離，γc為截止頻率。

1.3 濾波性能評(píng)價(jià)指標(biāo)

為了定量比較算法的性能，筆者采用等效視數(shù)（equivalent number of looks，ENL）作為本研究的評(píng)價(jià)標(biāo)準(zhǔn)，其計(jì)算方法為

式中：μ為圖像灰度均值，σ為圖像灰度標(biāo)準(zhǔn)方差。μ和σ分別用于描述算法對(duì)均值的保持能力和對(duì)方差的抑制能力。它們的比值ENL表現(xiàn)了圖形灰度的對(duì)比度，ENL越大，表示受到的干擾越小，整體的圖像效果越好。另外，本研究也考慮了算法在機(jī)器上對(duì)圖像處理速度的時(shí)間耗費(fèi)（Time cost）。

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

為驗(yàn)證本研究提出的模型可在保存紋理細(xì)節(jié)的同時(shí)盡可能地去除圖像中的斑點(diǎn)噪聲，減少達(dá)到穩(wěn)定時(shí)的迭代次數(shù)，實(shí)驗(yàn)選取了含有大量斑點(diǎn)噪聲的眼科高頻超聲圖像。

采用真實(shí)的UBM圖像，在實(shí)驗(yàn)中SRAD與改進(jìn)后的八方向SRAD時(shí)間步長(zhǎng)為0.05，迭代次數(shù)均為5次。在Win7系統(tǒng)下，采用Matlab2010b對(duì)圖像進(jìn)行算法的驗(yàn)證處理，結(jié)果如圖2所示。圖2A～D分別為原始噪聲圖像、直接采用SRAD濾波方式的圖像、采用CEDIF直接濾波的圖像、采用本研究方法濾波后重建的圖像。

從視覺角度來(lái)看，從圖2中可明顯看出，直接采用SRAD濾波可一定程度上抑制斑點(diǎn)噪聲，對(duì)圖像進(jìn)行平滑；采用CEDIF濾波出現(xiàn)斑塊效應(yīng)，影響圖像的真實(shí)性；而采用分層濾波后的圖像，斑點(diǎn)噪聲得到有效抑制，且在相同的迭代次數(shù)下，斑點(diǎn)噪聲抑制效果優(yōu)于原SRAD算法，圖像邊緣以及細(xì)節(jié)部分也得到了比較好的保持。

圖2 SRAD和CEDIF及本研究算法對(duì)比圖

表1 SRAD、CEDIF及Laplacian金字塔的多尺度各向異性斑點(diǎn)去噪算法的定量評(píng)估（±s）

表1 SRAD、CEDIF及Laplacian金字塔的多尺度各向異性斑點(diǎn)去噪算法的定量評(píng)估（±s）

注：“—”表示無(wú)此項(xiàng)；SRAD—各向異性斑點(diǎn)去噪；CEDIF—相干增強(qiáng)濾波

圖像原始圖像SRAD算法CEDIF算法Laplacian金字塔的多尺度各向異性斑點(diǎn)去噪均值50.305 9±0.581 4 49.621 8±0.494 4 45.352 9±0.698 3 50.521 8±0.285 4標(biāo)準(zhǔn)方差45.780 7±0.203 4 44.958 1±0.625 7 48.002 0±0.306 0 43.694 7±0.384 8等效視數(shù)1.096 0±0.034 9 1.086 2±0.042 2 0.950 6±0.029 7 1.156 3±0.025 9時(shí)間耗費(fèi)—1.397 0±0.015 0 1.475 7±0.015 8 0.945 9±0.017 0

為了量化比較，對(duì)相鄰5幀圖像進(jìn)行測(cè)量處理（表1）。由此可看出，本研究算法對(duì)圖像灰度均值的保持以及對(duì)方差的抑制能力優(yōu)于另外兩種算法，等效視數(shù)高于其他兩種算法，且在現(xiàn)有相同機(jī)器環(huán)境下，本研究算法的時(shí)間耗費(fèi)也比其他兩種小，且效率高。

3 結(jié)論

筆者提出的八方向SRAD-拉普拉斯金字塔模型能有效地抑制斑點(diǎn)噪聲，通過(guò)拉普拉斯金字塔，分離出子帶圖像；然后分析每層圖像不同的特點(diǎn)，再對(duì)圖像進(jìn)行側(cè)重點(diǎn)不同的處理，選取相應(yīng)子帶濾波方法，在有效抑制斑點(diǎn)噪聲的同時(shí)，對(duì)邊緣及細(xì)節(jié)進(jìn)行了很好的保護(hù)。通過(guò)對(duì)真實(shí)超聲圖像的分析證實(shí)，本研究提出的算法在保持邊緣紋理細(xì)節(jié)的同時(shí)，去噪效果也得到了加強(qiáng)，是一種快速有效的超聲圖像斑點(diǎn)抑制方法；另外該算法的空間復(fù)雜度對(duì)于其他算法而言并未太過(guò)復(fù)雜，無(wú)論軟件還是硬件上開銷均影響甚微，能夠迅速提升濾波速度，實(shí)用性較高[12]，將來(lái)對(duì)于便攜式超聲儀的研發(fā)也有一定的參考意義。

[1]楊軍,毛志林,何平,等.專用超聲診斷儀[M]//伍于添.醫(yī)學(xué)超聲設(shè)備原理設(shè)計(jì)應(yīng)用.北京:科學(xué)技術(shù)文獻(xiàn)出版社,2012. Yang J,Mao ZL,He P,et al.Special ultrasonic diagnostic instruments[M]//Wu YT.Medcine ultrasonic instruments princeple design application.Beijing:Scientific and Technical Documentation Press,2012.

[2]周盛,楊軍.全數(shù)字B超的信號(hào)處理[J].國(guó)際生物醫(yī)學(xué)工程雜志, 2007,30(6):356-360.DOI:10.3760/cma.j.issn.1673-4181.2007. 06.010. Zhou S,Yang J.Signal processing of totally digitalized B-ultrasonography[J].Int J Biomed Eng,2007,30(6):356-360.DOI: 10.3760/cma.j.issn.1673-4181.2007.06.010.

[3]Lee JS.Speckle suppression and analysis for synthetic aperture radar images[J].Opt Eng,1986,25(5):636-643.DOI:10.1117/12. 7973877.

[4]LimJS.Two-dimensionalsignalandimageprocessing[M]. Englewood Cliffs:Prentice-Hall,1990.

[5]Entrekin R,Jackson P,Jago JR,et al.Real time spatial compound imaginginbreastultrasound:technologyandearlyclinical experience[J].Medicamundi,1999,43(3):35-43.DOI:10.1016/ S0301-5629(01)00490-2.

[6]Burt PJ,Adelson EH.The Laplacian pyramid as a compact image code[J].IEEE T Commun,1983,31(4):532-540.DOI:10.1109/ TCOM.1983.1095851.

[7]楊先鳳,游書濤,彭博,等.醫(yī)學(xué)超聲圖像SRAD模型優(yōu)化研究[J].計(jì)算機(jī)仿真,2011,28(7):290-292.DOI:10.3969/j.issn.1006-9348. 2011.07.072. Yang XF,You ST,Peng B,et al.Optimization of SRAD model for medical ultrasound images[J].Computer Simulation,2011,28(7): 290-292.DOI:10.3969/j.issn.1006-9348.2011.07.072.

[8]Yu Y,Acton ST.Speckle reducing anisotropic diffusion[J].IEEE Trans Image Process,2002,11(11):1260-1270.DOI:10.1109/TIP. 2002.804276.

[9]Perona P,Malik J.Scale-space and edge detection using anisotropic diffusion[J].IEEE T Pattern Anal,1990,12(7):629-639.DOI:10. 1109/34.56205.

[10]Weickert J.Coherence-enhancing diffusion of colour images[J]. Image Vision Comput,1999,17(3/4):201-212.DOI:10.1016/S0262-8856(98)00102-4.

[11]Oppenheim AV,Schafer RW,Stockham TG.Nonlinear filtering of multiplied and convolved signals[J].IEEE Trans Audio and Electro, 1968,56(8):1264-1291.DOI:10.1109/PROC.1968.6570.

[12]王延群.走向國(guó)際市場(chǎng)的眼科超聲診療設(shè)備——記國(guó)產(chǎn)眼科超聲診療設(shè)備研發(fā)[J].國(guó)際生物醫(yī)學(xué)工程雜志,2012,35(1):1-2. DOI:10.3760/cma.j.issn.1673-4181.2012.01.001. Wang YQ.Ultrasonic ophthalmic medical equipment go out into international market—study on the research and development of the domestic ultrasonic diagnostic equipment[J].Int J Biomed Eng, 2012,35(1):1-2.DOI:10.3760/cma.j.issn.1673-4181.2012.01.001.

Research on speckle reduction in ophthalmology high frequency ultrasound imaging

Tian Ling,Zhou Sheng, Wang Xiaochun,Ye Qingsheng,Wang Yanqun
Institute of Biomedical Engineering,Chinese Academy of Medical Sciences&Peking Union Medical College,Tianjin 300192,China
Corresponding author:Wang Yanqun,Email:wyq@meda.com.cn

ObjectiveSpeckle is the inherent problem exists in B-Mode ultrasound image,and effective speckle noise reduction will improve the image quality and contribute greatly to clinical doctors to diagnose, especially in fine ophthalmic examination with high frequency ultrasound.MethodsThis paper proposed a new speckle reduction method based on the Laplacian pyramid and multiscale analysis.In the Laplacian pyramid,true clinical features were separated from noise,according to the different bandpass ultrasound image characteristics in each layer,and the advanced eight directions speckle reduction anisotropic diffusion(SRAD)was adapted to suppressed the speckle noise,and the identified noise and the extracted features were selectively emphasized by suitable edge,coherence and contrast enhancement filtering from fine to coarse scales.The performance of the proposed method was compared to the traditional SRAD method and coherence-enhancing diffusion method by measuring the equivalent number of looks(ENL)and Time cost.ResultsThe ENL and Time cost value of the proposed method were higher compared to the SRAD and the cedif method,i.e 1.172 3 vs 1.122 3,0.929 3 and 0.864 0 vs 1.396 0,1.468 3.ConclusionsIn summary,the proposed method can more effectively remove the speckle noise while preserving the edge and details of the images.

Speckle noise;Laplacian pyramid;Multiscale analysis;Anisotropic diffusion

王延群，Email:wyq@meda.com.cn

10．3760/cma．j．issn．1673-4181．2016．01．002

國(guó)家科技支撐計(jì)劃課題（2012BAI13B02）；中央級(jí)公益性科研院所基本科研業(yè)務(wù)專項(xiàng)課題

2015-11-30）