二進(jìn)小波與擴(kuò)散濾波結(jié)合的光學(xué)相干層析圖像降噪

張 田，孫延奎，田小林

(1.清華大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)系，北京100084;2.澳門科技大學(xué)資訊科技學(xué)院，澳門特別行政區(qū))

其中q0為待處理圖像中的斑點(diǎn)方差系數(shù)，ρ是常量，t=n×Δt，n是迭代次數(shù)，Δt是迭代步長。

與離散小波變換相比，二進(jìn)小波變換具有平移不變性，因此，二進(jìn)小波在圖像降噪中具有重要應(yīng)用［10］。二進(jìn)小波變換具有尺度相關(guān)性，即圖像重要邊緣的小波系數(shù)會(huì)分布在多個(gè)尺度中，而圖像噪聲的小波系數(shù)主要集中在最低尺度中;特別的，圖像重要邊緣處的小波系數(shù)幅值隨著尺度升高不會(huì)明顯減小;而圖像噪聲的小波系數(shù)幅值隨著尺度升高會(huì)急劇衰減。

光學(xué)相干層析成像技術(shù)［1］(OCT:Optical Coherence Tomography)是一種基于弱相干光的干涉原理，通過檢測不同組織對入射的弱相干光的背向反射或散射信號得到生物組織的二維圖像或三維結(jié)構(gòu)的光學(xué)成像技術(shù)。因?yàn)槠渚哂谐上窨?、分辨率高、非侵入式等?yōu)點(diǎn)，自從上世紀(jì)90年代以來，OCT在醫(yī)學(xué)領(lǐng)域獲得了廣泛應(yīng)用，尤其在眼科疾病診斷方面，OCT系統(tǒng)獲取的高清晰(μm量級)視網(wǎng)膜圖像對于青光眼等眼科疾病的診斷提供了極大的幫助［2］。

但在OCT成像設(shè)備的實(shí)際應(yīng)用中，光源、檢測電路、振鏡都會(huì)帶來噪聲，此外，光的多次散射也會(huì)產(chǎn)生很多散斑，這些都嚴(yán)重降低了圖像質(zhì)量，給OCT成像設(shè)備的臨床應(yīng)用帶來了不利影響。OCT圖像中存在的多種噪聲中，由于使用的光源的相干作用產(chǎn)生的斑點(diǎn)噪聲［3］是最主要的部分，這種噪聲以乘法噪聲的形式存在，其數(shù)學(xué)模型

乘法噪聲的去除對降噪方法提出了新的要求，多年來國內(nèi)外研究者對這個(gè)問題進(jìn)行了大量研究，取得了眾多成果。一類是對原圖像進(jìn)行一次對數(shù)運(yùn)算，將乘法噪聲轉(zhuǎn)化為加法噪聲后，再應(yīng)用如BiShrink［4］，ProbShrink［5］，曲波(curvelet)系數(shù)收縮方法［6］等傳統(tǒng)加法噪聲去除方法。另一類是利用斑點(diǎn)噪聲的特性開發(fā)的專用于去除乘法斑點(diǎn)噪聲的方法，如利用歸一化小波系數(shù)的小波擴(kuò)散降斑方法［7］，斑點(diǎn)抑制各向異性擴(kuò)散(SRAD: Speckle Reduction Anisotropic Diffusion)［8］等。SRAD算法是在傳統(tǒng)各向異性擴(kuò)散濾波方法［9］的基礎(chǔ)上，通過設(shè)計(jì)一種專用于檢測乘法斑點(diǎn)噪聲的邊緣檢測算子，取得了較好的抑制斑點(diǎn)噪聲的效果。但在對視網(wǎng)膜OCT圖像的降噪實(shí)驗(yàn)中發(fā)現(xiàn)，SRAD算法雖然具有較強(qiáng)的去除斑點(diǎn)噪聲的能力，但仍有很多較強(qiáng)的噪聲會(huì)被誤判為邊緣從而被保留，使得SRAD降噪后的圖像中仍存在較大量噪聲。而小波變換具有很強(qiáng)的分離噪聲和邊緣的能力，而且圖像的噪聲和邊緣的小波系數(shù)幅值在不同尺度間具有不同的衰減性。因此通過結(jié)合SRAD的降噪原理和小波系數(shù)尺度間的相關(guān)性，筆者提出一種改進(jìn)的SRAD算法，在保持SRAD增強(qiáng)邊緣的能力的同時(shí)，進(jìn)一步提高其降噪能力。

1 算法描述

筆者利用OCT圖像的對數(shù)圖像的二進(jìn)小波系數(shù)幅值的尺度相關(guān)性對邊緣還是均質(zhì)區(qū)域的指示作用，改進(jìn)SRAD算法中的擴(kuò)散系數(shù)計(jì)算公式，從而增強(qiáng)其保持邊緣并去除斑點(diǎn)的能力。因?yàn)樾路椒ńY(jié)合了二進(jìn)小波和SRAD算法的優(yōu)點(diǎn)，故稱為 DW－SRAD(Dyadic Wavelet Guided Speckle Reducing Anisotropic Diffusion)算法。

1.1 斑點(diǎn)抑制擴(kuò)散濾波算法［8］

SRAD算法是在各向異性擴(kuò)散濾波算法的基礎(chǔ)上針對超聲圖像中普遍存在的乘性斑點(diǎn)噪聲開發(fā)的。其具體過程如下:

對于一張能量有限的帶有斑點(diǎn)噪聲的圖像I0(x，y)，令其支撐區(qū)間為Ω，則可根據(jù)下列方程計(jì)算其每一步迭代擴(kuò)散后得到的新圖像:

其中?Ω代表支撐區(qū)間Ω的邊界，n為?Ω的外法向量，c(q)是擴(kuò)散系數(shù)，

其中q(x，y;t)為方差系數(shù)。

q0(t)為斑點(diǎn)度量函數(shù)。

其中q0為待處理圖像中的斑點(diǎn)方差系數(shù)，ρ是常量，t=n×Δt，n是迭代次數(shù)，Δt是迭代步長。

1.2 二進(jìn)小波變換及其尺度相關(guān)性

與離散小波變換相比，二進(jìn)小波變換具有平移不變性，因此，二進(jìn)小波在圖像降噪中具有重要應(yīng)用［10］。二進(jìn)小波變換具有尺度相關(guān)性，即圖像重要邊緣的小波系數(shù)會(huì)分布在多個(gè)尺度中，而圖像噪聲的小波系數(shù)主要集中在最低尺度中;特別的，圖像重要邊緣處的小波系數(shù)幅值隨著尺度升高不會(huì)明顯減小;而圖像噪聲的小波系數(shù)幅值隨著尺度升高會(huì)急劇衰減。

圖像I(x，y)連續(xù)做2次二維二進(jìn)小波變換后，可得到一個(gè)低頻圖像和4個(gè)高頻圖像(x，y)(x，y)(x，y)，(x，y)。在二維小波域，每一個(gè)點(diǎn)(x0，y0)在第j尺度都有一個(gè)二維影響錐(COI:Cone of Influence(x0，y0)［10－12］。令N2j/(x0，y0)表示影響錐(x0，y0)中各小波系數(shù)絕對值之和，即

定義:

稱τ(x0，y0)為點(diǎn)(x0，y0)處二進(jìn)小波系數(shù)尺度間的相關(guān)系數(shù)。二進(jìn)小波變換的尺度相關(guān)性表明，較大的τ(x0，y0)意味著點(diǎn)(x0，y0)更有可能位于重要邊緣上;較小的τ(x0，y0)意味著點(diǎn)(x0，y0)很可能位于均質(zhì)區(qū)域。這表明，τ(x0，y0)能夠起到指示點(diǎn)(x0，y0)是位于圖像的重要邊緣上還是位于噪聲的均質(zhì)(homogenous)區(qū)域的作用。利用τ值的特征，可改進(jìn)SRAD中的擴(kuò)散系數(shù)計(jì)算公式，從而進(jìn)一步提高SRAD降低斑點(diǎn)噪聲并保持OCT圖像邊緣的能力。

1.3 算法流程

筆者提出了一個(gè)結(jié)合各向異性擴(kuò)散濾波算法和二進(jìn)小波系數(shù)尺度間相關(guān)性的OCT圖像降斑降噪方法。該方法利用OCT圖像的對數(shù)圖像在二進(jìn)小波域的小波系數(shù)尺度間的相關(guān)系數(shù)，設(shè)計(jì)一種更適合OCT圖像的擴(kuò)散系數(shù)計(jì)算公式，從而起到更好的保留OCT圖像重要邊緣，同時(shí)有效去除斑點(diǎn)噪聲的效果。具體操作過程如下:

1)計(jì)算二進(jìn)小波系數(shù)尺度間的相關(guān)性系數(shù)。

對原圖像I做一次對數(shù)運(yùn)算，得到對數(shù)圖像Ilog;對Ilog進(jìn)行2次二維二進(jìn)小波變換，得到一個(gè)低頻圖像S22(x，y)和4個(gè)高頻圖像(x，y)，(x，y)，(x，y)，(x，y)。再根據(jù)式(7)、式(8)計(jì)算每點(diǎn)處二進(jìn)小波系數(shù)尺度間的相關(guān)性系數(shù)τ(x0，y0)。

2)設(shè)計(jì)新的擴(kuò)散系數(shù)計(jì)算公式。

利用τ具有區(qū)分噪聲和邊緣的能力，引入一個(gè)擴(kuò)散系數(shù)修正因子s(τ):定義s(τ)為τ的單調(diào)遞減函數(shù)，使得在圖像的邊緣上，τ相對較大時(shí)s(τ)較小;而在圖像的均質(zhì)區(qū)域，τ較小時(shí)s(τ)較大。例如，可定義s(τ)的計(jì)算公式為:

其中K為一個(gè)常量，由實(shí)驗(yàn)給出。

由s(τ)修正SRAD的擴(kuò)散系數(shù)計(jì)算公式c(q)得到一個(gè)新的擴(kuò)散系數(shù)計(jì)算公式。

3)利用cnew(q，τ)代替SRAD中的擴(kuò)散系數(shù)計(jì)算公式c(q)進(jìn)行擴(kuò)散濾波

與c(q)相比，cnew(q，τ)具有更好的定位邊緣與噪聲的能力，因而將cnew(q，τ)應(yīng)用于圖像的迭代擴(kuò)散降噪過程，會(huì)具有更好的降斑降噪能力，即在避免對邊緣區(qū)域過平滑的同時(shí)更有效地抑制斑點(diǎn)噪聲。

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

為驗(yàn)證所提出算法的有效性，筆者利用視網(wǎng)膜OCT圖像進(jìn)行實(shí)驗(yàn)。數(shù)據(jù)采集系統(tǒng)的有關(guān)參數(shù)為:光源中心波長 830 nm，有效探測深度3.4 mm，信噪比為51 dB，縱向分辨率8.5μm［13］。實(shí)驗(yàn)中，迭代步長Δt為0.25，ρ取1，迭代次數(shù)為50次，常數(shù)K為20。對15張視網(wǎng)膜OCT圖像進(jìn)行了降噪實(shí)驗(yàn)，其中一幅圖像及其降噪后圖像的對比結(jié)果如圖1所示。

圖1 降噪前后的效果對比圖Fig.1 The denoising results of a retina OCT image by DW_SRAD method

圖1(a)為原OCT圖像，綠色框中為背景區(qū)域，紅色框中為前景區(qū)域，用于計(jì)算4種降噪指標(biāo);圖1(b)為圖1(a)中藍(lán)線標(biāo)注的A－scan中的信號;圖1(c)為DW_SRAD降噪后的圖像;圖1(d)為圖1(c)中藍(lán)線標(biāo)注的A－scan中的信號。從圖1中可看出，原圖中大量噪聲被有效的降低了，包括黑色背景區(qū)域中的離散噪聲和高亮前景區(qū)域中的噪聲。圖1中(b)和(d)的對比能清楚地看到重要邊緣信號得到了較好的保留。

為定量地評價(jià)算法降噪效果，該次實(shí)驗(yàn)引入了幾個(gè)國際通行的指標(biāo):

1)信噪比SNR(Signal－to－Noise Ratio)［14］

用于衡量圖像中信號和噪聲強(qiáng)度:

其中μ為信號強(qiáng)度的期望值;σ為噪聲的標(biāo)準(zhǔn)差。

2)等效視數(shù) ENL(Effective Number of Looks)［14］

用于度量圖像中由于噪聲引入的擾動(dòng)，即一塊本應(yīng)該是平滑區(qū)域但被噪聲污染后的區(qū)域平滑性度量

其中t代表圖像中的前景區(qū)域，即圖1中紅色框中的區(qū)域，μt為該區(qū)域中的像素均值，σt為該區(qū)域中像素標(biāo)準(zhǔn)差。

3)對比度噪聲比率 CNR(Contrast－to－Noise Ratio)［14］

提供圖像特征區(qū)域(感興趣區(qū)域)與背景區(qū)域(例如空氣區(qū)域)的有效對比度的一個(gè)客觀度量，計(jì)算公式:

其中符號t，μt和σt的含義與式(12)中的相同。r表示圖像中的背景參考區(qū)域，即圖1中的綠色框中的區(qū)域，μr和σr分別為該區(qū)域中的信號均值與標(biāo)準(zhǔn)差。

4)邊緣保持參數(shù)β［15］

用于度量降噪后的圖像邊緣銳度降低的程度。

其中I和Id分別為原OCT圖像和降噪后的OCT圖像，Γ(X，Y)=∑XijYij，ΔI和ΔId分別為對原圖像和降噪后圖像進(jìn)行拉普拉斯濾波的結(jié)果和為圖像像素均值，該實(shí)驗(yàn)中使用3×3拉氏算子。

對15張視網(wǎng)膜OCT圖像利用SARD與DW_ SRAD兩種算法進(jìn)行實(shí)驗(yàn)后，得到的以上4種指標(biāo)的平均結(jié)果如表1所示。

表1 DW_SRAD與SRAD降噪性能對比Table 1 Performance indexes comparison between the DW_SRAD and SRAD method

從表1可看出，新方法在沒有明顯加重對原圖像的模糊(用于衡量降噪后圖像銳度或者與原圖像相比的失真程度的指標(biāo)β值保持基本不變)的前提下，其它3種衡量降噪效果的指標(biāo)都有明顯的提高。

3 結(jié)論

利用相鄰尺度間的二進(jìn)小波系數(shù)的相關(guān)性，提出了一種新的擴(kuò)散系數(shù)計(jì)算公式，對斑點(diǎn)抑制各向異性擴(kuò)散濾波方法SRAD進(jìn)行了改進(jìn)。改進(jìn)后的算法在保持經(jīng)典SRAD算法的增強(qiáng)邊緣能力的同時(shí)，顯著提高了其降低噪聲的能力，信噪比、等效視覺等客觀量化指標(biāo)都有所提高。利用小波系數(shù)的尺度相關(guān)性可以改善擴(kuò)散濾波效果。該方法對一般斑點(diǎn)圖像噪聲降噪有普適性。

［1］Zysk A M，Nguyen F T，Oldenburg A L，et al.Optical coherence tomography:a review of clinical development from bench to bedside［J］.Journal of Biomedical Optics，2007，12(5):12－16.

［2］Jindahra P，Hedges T R，Mendoza－Santiesteban C E，et al.Optical coherence tomography of the retina:applications in neurology［J］.Current Opinion In Neurology，2010，23(1):16－23.

［3］Karamata B，Hassler K，Laubscher M，etal.Speckle statistics in optical coherence tomography［J］.Journal of the Optical Society of America A－Optics Image Science and Vision，2005，22(4):593－596.

［4］Sendur L，Selesnick IW.Bivariate shrinkage functions for wavelet－based denoising exploiting interscale dependency［J］.IEEE Transactions on Signal Processing，2002，50(11):2744－2756.

［5］Pizurica A，Philips W.Estimating the probability of the presence of a signal of interest inmultiresolution singleand multiband image denoising［J］.IEEE Transactions on Image Processing，2006，15(3):654－665.

［6］Jian Zhong－ping，Yu Zhao－xia，Yu Ling－feng，et al. Speckle attenuation in optical coherence tomography by curvelet shrinkage［J］.Optics Letters，2009，34(10): 1516－1518.

［7］Yong Yue，Croitoru M M，Bidani A，et al.Nonlinear multiscale wavelet diffusion for speckle suppression and edge enhancement in ultrasound images［J］.IEEE Transactions on Medical Imaging，2006，25(3):297－311.

［8］Yu Y J，Acton S T.Speckle reducing anisotropic diffu－sion［J］.IEEE Transactions on Image Processing，2002，11(11):1260－1270.

［9］Perona P，Malik J.Scale－Space and edge－detection using anisotropic diffusion［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1990，12(7):629－639.

［10］Mallat S，Zhong S.Singularity detection and processing with wavelets［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1992，14(7):710－732.

［11］Junmei Zhong，Ruola Ning，David Conover.Image denoising based on multiscale singularity detection for cone beam CT breast imaging［J］.IEEE Transactions on Medical Imaging，2004，23(6):696－703.

［12］Zhong Jun－mei，Sun Hui－fang.Wavelet－based multiscale anisotropic diffusion with Adaptive Statistical analysis for image restoration［J］.IEEE Transactions on Circuits and Systems—I:Regular Papers，2008，55(9):2716－2725.

［13］段煉，何永紅，朱銳，等.三維譜域光學(xué)相干層析成像系統(tǒng)的研制［J］.中國激光，2009，36(10):2528－2533.

Duan Lian，He Yong－h(huán)ong，Zhu Rui，et al.Development of a spectrum domain 3D optical coherence tomography system［J］.Chinese Journal of Lasers，2009，36(10): 2528－2533.

［14］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 Transactions on Medical Imaging，2007，26(6): 761－771.

［15］鄧菊香，梁艷梅.光學(xué)相干層析圖像的小波降噪方法研究［J］.光學(xué)學(xué)報(bào)，2009，29(8):2138－2141.

Deng Ju－xiang，Liang Yan－mei.Noise reduction with wavelet transform in optical coherence tomographic images［J］.Acta Optica Sinica，2009，29(8):2138－2141.