二進小波與擴散濾波結合的光學相干層析圖像降噪

張 田，孫延奎，田小林

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

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

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

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

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

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

1 算法描述

筆者利用OCT圖像的對數(shù)圖像的二進小波系數(shù)幅值的尺度相關性對邊緣還是均質(zhì)區(qū)域的指示作用，改進SRAD算法中的擴散系數(shù)計算公式，從而增強其保持邊緣并去除斑點的能力。因為新方法結合了二進小波和SRAD算法的優(yōu)點，故稱為 DW－SRAD(Dyadic Wavelet Guided Speckle Reducing Anisotropic Diffusion)算法。

1.1 斑點抑制擴散濾波算法［8］

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

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

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

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

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

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

1.2 二進小波變換及其尺度相關性

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

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

定義:

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

1.3 算法流程

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

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

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

2)設計新的擴散系數(shù)計算公式。

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

其中K為一個常量，由實驗給出。

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

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

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

2 實驗結果

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

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

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

為定量地評價算法降噪效果，該次實驗引入了幾個國際通行的指標:

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

用于衡量圖像中信號和噪聲強度:

其中μ為信號強度的期望值;σ為噪聲的標準差。

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

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

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

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

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

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

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

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

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

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

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

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

3 結論

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

［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］段煉，何永紅，朱銳，等.三維譜域光學相干層析成像系統(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］鄧菊香，梁艷梅.光學相干層析圖像的小波降噪方法研究［J］.光學學報，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.