張 田，孫延奎，田小林

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

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

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

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

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

1 算法描述

筆者利用OCT圖像的對(duì)數(shù)圖像的二進(jìn)小波系數(shù)幅值的尺度相關(guān)性對(duì)邊緣還是均質(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ǔ)上針對(duì)超聲圖像中普遍存在的乘性斑點(diǎn)噪聲開(kāi)發(fā)的。其具體過(guò)程如下:

對(duì)于一張能量有限的帶有斑點(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是迭代步長(zhǎng)。

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ù)絕對(duì)值之和，即

定義:

稱τ(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圖像的對(duì)數(shù)圖像在二進(jìn)小波域的小波系數(shù)尺度間的相關(guān)系數(shù)，設(shè)計(jì)一種更適合OCT圖像的擴(kuò)散系數(shù)計(jì)算公式，從而起到更好的保留OCT圖像重要邊緣，同時(shí)有效去除斑點(diǎn)噪聲的效果。具體操作過(guò)程如下:

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

對(duì)原圖像I做一次對(duì)數(shù)運(yùn)算，得到對(duì)數(shù)圖像Ilog;對(duì)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ù)，使得在圖像的邊緣上，τ相對(duì)較大時(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ò)散降噪過(guò)程，會(huì)具有更好的降斑降噪能力，即在避免對(duì)邊緣區(qū)域過(guò)平滑的同時(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ù)為:光源中心波長(zhǎng) 830 nm，有效探測(cè)深度3.4 mm，信噪比為51 dB，縱向分辨率8.5μm［13］。實(shí)驗(yàn)中，迭代步長(zhǎng)Δt為0.25，ρ取1，迭代次數(shù)為50次，常數(shù)K為20。對(duì)15張視網(wǎng)膜OCT圖像進(jìn)行了降噪實(shí)驗(yàn)，其中一幅圖像及其降噪后圖像的對(duì)比結(jié)果如圖1所示。

圖1 降噪前后的效果對(duì)比圖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中的信號(hào);圖1(c)為DW_SRAD降噪后的圖像;圖1(d)為圖1(c)中藍(lán)線標(biāo)注的A－scan中的信號(hào)。從圖1中可看出，原圖中大量噪聲被有效的降低了，包括黑色背景區(qū)域中的離散噪聲和高亮前景區(qū)域中的噪聲。圖1中(b)和(d)的對(duì)比能清楚地看到重要邊緣信號(hào)得到了較好的保留。

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

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

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

其中μ為信號(hào)強(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)對(duì)比度噪聲比率 CNR(Contrast－to－Noise Ratio)［14］

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

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

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

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

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

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

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

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

3 結(jié)論

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