基于貝葉斯最小均方差的OCT圖像散斑處理

, , ,

(上海師范大學(xué) 信息與機(jī)電工程學(xué)院,上海 200234)

光學(xué)相干層析技術(shù)(OCT)是近20年來(lái)發(fā)展快速的一種生物醫(yī)學(xué)層析成像技術(shù).由于OCT是基于相干光探測(cè)的成像技術(shù),在圖像中不可避免的會(huì)產(chǎn)生散斑.這種散斑噪聲降低了圖像的質(zhì)量,掩蓋了生物組織的結(jié)構(gòu)細(xì)節(jié),給醫(yī)學(xué)診斷帶來(lái)了很大的困難,所以對(duì)OCT圖像散斑噪聲處理有著很大的現(xiàn)實(shí)意義.

現(xiàn)在國(guó)內(nèi)外對(duì)散斑處理的研究主要基于兩個(gè)方面.一方面是改進(jìn)OCT成像系統(tǒng).OCT的成像質(zhì)量與光源、數(shù)據(jù)采集速度、掃描速度、測(cè)量靈敏度等性能參數(shù)有關(guān),Kim[1]等人提出了利用部分空間相干寬帶光源來(lái)抑制散斑噪聲.Pircher[2]等人提出了頻域合成法,利用兩個(gè)非相干信號(hào)來(lái)減少散斑,同時(shí)還保留了較高的空間分辨率.Desjardins[3]等人通過(guò)在多重散射角度上使用角度合成,能夠降低散斑噪聲等等.然而,以上這些方法都需要重新設(shè)計(jì)OCT系統(tǒng),數(shù)據(jù)采集的時(shí)間也較長(zhǎng),并且在多數(shù)情況下散斑的抑制效果還受到一定程度的限制.

另一方面是研究OCT圖像降噪算法.散斑抑制的算法仍然是OCT圖像后處理的一個(gè)重要部分.通過(guò)對(duì)散斑噪聲頻譜分布、統(tǒng)計(jì)特性、圖像特點(diǎn)的分析,提出了很多抑制OCT圖像散斑的方法.Pizurica[4]等人提出了基于小波變換的方法,該方法減少了散斑噪聲,但降低了微小結(jié)構(gòu)細(xì)節(jié)的可見(jiàn)度;Puvanathasan和Bizheva[5]提出的Type II Fuzzy AD方法對(duì)低噪聲的去除有好的效果,但對(duì)高噪聲的濾除有一定的局限性.柯麗[6]等人提出了應(yīng)用多級(jí)維納濾波的方法對(duì)OCT圖像去噪,該方法較好的抑制了散斑噪聲,并且圖像的視覺(jué)效果良好.蘇鵬[7]等人提出了雙樹(shù)復(fù)小波和混合概率模型的方法抑制OCT圖像的散斑,該方法抑制了散斑噪聲并保持圖像邊緣銳度的相對(duì)穩(wěn)定性.

1 OCT圖像中的散斑

1.1 散斑的來(lái)源

散斑是由隨機(jī)相位的散射光波干涉疊加產(chǎn)生的.在OCT成像系統(tǒng)中,光照射到存在大量的散射顆粒的樣品組織內(nèi)部發(fā)生散射,只有與參考光束光程差在相干長(zhǎng)度范圍內(nèi)的散射光才能與參考光相干成像,然而在高散射生物組織中,相干長(zhǎng)度內(nèi)存在大量的散射顆粒截面,而被光電探測(cè)器接收的干涉光中既有單次背向散射光,又有多次散射光,光電探測(cè)器上會(huì)有光程差為nπλ的相干散射光束同時(shí)到達(dá),產(chǎn)生高斯包絡(luò)的具有相位差為nπ的交變電信號(hào),它們彼此相干疊加就形成了散斑.不僅大量散射顆粒形成多次散射,在相干長(zhǎng)度內(nèi)不同深度截面上的散射光由于光程差引起的相位差也會(huì)形成散斑.另外生物組織還會(huì)吸收一部分光,使光的一部分頻率丟失,這些頻率的變化也會(huì)導(dǎo)致散斑的出現(xiàn).所以,散斑的形成是非常復(fù)雜而且不可避免的[8-11].

1.2 散斑的數(shù)學(xué)模型

在OCT層析圖像中,散斑依賴于成像光束的波長(zhǎng)和成像對(duì)象的結(jié)構(gòu)細(xì)節(jié).因此,散斑攜帶了兩個(gè)信息:成像對(duì)象的結(jié)構(gòu)和噪聲分量.由散斑噪聲的統(tǒng)計(jì)特性可知,散斑噪聲是一個(gè)乘性噪聲[12],因此 OCT圖像的信號(hào)模型可以表示為:

s(x,y)=f(x,y)·n(x,y) .

(1)

其中,(x,y)是OCT圖像上的一個(gè)像素點(diǎn)二維坐標(biāo),s(x,y)是檢測(cè)到的OCT圖像信號(hào)數(shù)據(jù),f(x,y)表示理想的無(wú)噪聲信號(hào)數(shù)據(jù),n(x,y)表示散斑噪聲信號(hào)數(shù)據(jù).

設(shè)M表示OCT圖像的像素點(diǎn)的集合,m∈M表示OCT圖像上的一個(gè)像素點(diǎn)的二維坐標(biāo)(x,y),則表達(dá)式(1)可以表示成如下形式:

s(m)=f(m)·n(m) .

(2)

由于散斑噪聲是一個(gè)乘性噪聲,從散斑噪聲中分離出無(wú)噪聲信號(hào)的問(wèn)題非常困難.可將檢測(cè)到的OCT圖像信號(hào)數(shù)據(jù)放在對(duì)數(shù)空間中,使得乘性散斑噪聲變成一個(gè)加性噪聲,從而利于散斑噪聲的去除.將(2)式兩邊取對(duì)數(shù):

log[s(m)]=log[f(m)]+log[n(m)] .

(3)

即：

S(m)=F(m)+N(m) .

(4)

其中,

S(m)=log[s(m)]，F(xiàn)(m)=log[f(m)]，N(m)=log[n(m)] .

2 算法的基本思想

2.1 貝葉斯估計(jì)

本文作者提出了一種基于非線性對(duì)數(shù)空間的貝葉斯最小均方差估計(jì)[13]的算法.該算法的關(guān)鍵思想是根據(jù)高斯噪聲分布的特點(diǎn),從噪聲的高斯分布中抽取樣本,對(duì)樣本內(nèi)的像素點(diǎn)賦予相應(yīng)的權(quán)值,用加權(quán)直方圖來(lái)估計(jì)后驗(yàn)分布p(F(m)|S(m)),從而計(jì)算出無(wú)噪聲數(shù)據(jù)F(m)的貝葉斯估計(jì)值.

(5)

∑p(F(m)|S(m))F(m) .

(6)

表達(dá)式(6)顯示理想的無(wú)噪聲信號(hào)F(m)的最優(yōu)貝葉斯估計(jì),是以檢測(cè)到的OCT圖像信號(hào)的數(shù)據(jù)S(m)為條件的統(tǒng)計(jì)平均值.

2.2 圖像像素的權(quán)值

如果把圖像信號(hào)S看成是一個(gè)隨機(jī)變量,那么加性噪聲服從高斯分布.在對(duì)乘性散斑噪聲取對(duì)數(shù)之后,散斑噪聲轉(zhuǎn)變成了加性高斯噪聲分布.設(shè)m′是M中的一個(gè)隨機(jī)點(diǎn),m是圖像鄰域的中心像素點(diǎn),以像素點(diǎn)m為中心的高斯分布表示如下:

(7)

其中,‖m′-m‖表示像素點(diǎn)m′到中心像素點(diǎn)m的平方歐氏距離,σspatial是高斯分布的空間尺度因子,確定圖像鄰域的尺度大小.設(shè)置σspatial為7個(gè)像素,即圖像鄰域的模板大小為7×7.從高斯分布Q(m′|m)產(chǎn)生的點(diǎn)m′是以m點(diǎn)為中心的鄰域內(nèi)的點(diǎn).根據(jù)高斯分布的特點(diǎn),m′滿足的抽樣條件是:

|u(m)-u(m′)|<2σ.

(8)

一幅圖像鄰域內(nèi)的相鄰像素之間是高度相關(guān)的,某一點(diǎn)的灰度值與其周圍點(diǎn)的灰度值非常接近,在一幅圖像中,如果一個(gè)像素點(diǎn)的灰度值遠(yuǎn)大于或小于其鄰域的灰度值,則表明該像素點(diǎn)與其鄰域像素點(diǎn)的相關(guān)性小,該像素點(diǎn)很可能為噪聲點(diǎn).像素之間的相關(guān)程度可用高斯分布的似然函數(shù)來(lái)描述:

(9)

2.3 加權(quán)直方圖

因?yàn)楹篁?yàn)分布模型未知,所以這里采用直方圖估計(jì)后驗(yàn)分布.直方圖是對(duì)圖像或圖像的某個(gè)區(qū)域中每個(gè)灰度級(jí)分布情況的統(tǒng)計(jì).但是直方圖只提供了圖像的灰度分布信息,而圖像中的有效信號(hào)數(shù)據(jù)和噪聲數(shù)據(jù)的分布情況卻不知.為了區(qū)分圖像中的有效信號(hào)和噪聲,將圖像像素的權(quán)值引入到直方圖中,即加權(quán)直方圖.

加權(quán)直方圖[14-15]的思想是:在計(jì)算直方圖時(shí),給每個(gè)像素點(diǎn)賦予一定的權(quán)值,統(tǒng)計(jì)落入每一個(gè)直方圖區(qū)間像素的加權(quán)個(gè)數(shù).利用加權(quán)直方圖來(lái)計(jì)算后驗(yàn)分布p(F(m)|S(m)),p(F(m)|S(m))表示以檢測(cè)到的OCT圖像信號(hào)的數(shù)據(jù)S(m)為條件的無(wú)噪聲數(shù)據(jù)F(m)的概率分布.

(10)

其中,δ是 Delta函數(shù),z是歸一化常數(shù).根據(jù)歸一化條件:

可得:

圖1 算法流程圖

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

為了驗(yàn)證本算法的有效性和優(yōu)越性,本文作者使用OCT實(shí)時(shí)成像圖片進(jìn)行處理,分別運(yùn)用小波去噪,中值濾波去噪和本算法進(jìn)行降噪處理.通過(guò)實(shí)驗(yàn)數(shù)據(jù)比較,來(lái)展現(xiàn)本算法的優(yōu)越性.除了圖片的視覺(jué)對(duì)比外,作者選擇了2個(gè)定量的圖像質(zhì)量指標(biāo)作為客觀評(píng)價(jià)依據(jù),一個(gè)是信噪比,另一個(gè)是等效視數(shù).信噪比是用于比較被評(píng)價(jià)圖像與原圖像質(zhì)量的參數(shù),信噪比的數(shù)值越大,圖像中的散斑噪聲越少,圖像質(zhì)量越好.等效視數(shù)越大,表示圖像上的相干斑越弱,散斑對(duì)圖像的影響越小.信噪比[5](SAR)和等效視數(shù)[16](ENL)被定義如下:

SNR=10log10[max(S2)/σ2] .

(11)

(12)

其中,S代表圖像像素的灰度值,σ代表圖像像素的灰度值的標(biāo)準(zhǔn)差,u代表圖像像素灰度值的平均值.

選取2張膀胱OCT實(shí)時(shí)成像圖片,第一張圖片的實(shí)驗(yàn)結(jié)果對(duì)比圖像如圖2所示.其中,(a)為原始圖像,(b)為小波變換去噪后的結(jié)果圖,(c)為中值濾波去噪后的結(jié)果圖,(d)為本算法去噪后的結(jié)果圖,表1為其客觀對(duì)比數(shù)據(jù).

第二張圖片的實(shí)驗(yàn)結(jié)果對(duì)比圖像如圖3所示.其中,(1)為原始圖像,(2)為小波變換去噪后的結(jié)果圖,(3)為中值濾波去噪后的結(jié)果圖,(4)為本算法去噪后的結(jié)果圖,表2為其客觀對(duì)比數(shù)據(jù).

表1 OCT實(shí)時(shí)成像圖片的客觀數(shù)據(jù)

圖2 第一張OCT實(shí)時(shí)成像圖片處理結(jié)果

去噪方法SNRENL原始圖像11.78630.4598小波變換11.85720.4702中值濾波12.32940.4762本算法12.54350.4767

圖3 第二張OCT實(shí)時(shí)成像圖片處理結(jié)果

4 結(jié)束語(yǔ)

為了抑制OCT圖像的散斑噪聲,通過(guò)對(duì)OCT圖像散斑來(lái)源進(jìn)行分析,并根據(jù)OCT的成像特點(diǎn)以及散斑噪聲模型的分析,本文作者提出了一個(gè)基于對(duì)數(shù)空間的貝葉斯最小均方差估計(jì)的新算法.從表1和表2中的數(shù)據(jù)分析可知,與小波變換,中值濾波相比,算法在信噪比和等效視數(shù)方面都有明顯的改善.這表明本算法去除OCT圖像散斑后的圖像質(zhì)量更好,圖像上的相干斑更弱,散斑對(duì)圖像的影響更小.綜上所述,本算法有助于降低OCT圖像的散斑噪聲,提高OCT圖像的質(zhì)量,在醫(yī)學(xué)診斷中能夠?yàn)獒t(yī)護(hù)人員提供清晰準(zhǔn)確的圖像信息.

參考文獻(xiàn):

[1] KIM J,MILLER D,KIM E,et al.Optical coherence tomography speckle reduction by a partially spatially coherent source[J].Journal of Biomedical Optics,2005,10(6):1-9.

[2] PIRCHER M,GTZINGER E,LEITGEB R,et al.Speckle reduction in optical coherence tomography by frequency compounding [J].Journal of Biomedical Optics,2003,8(3):565-569.

[3] DESJARDINS A E,VAKOC B J,OH W Y.Angle-resolved Optical Coherence Tomography with sequential angular selectivity for speckle reduction[J].Optics Express,2007,15(10):6200-6209 .

[4] PIZURICA A,PHILIPS W,LEMAHIEU I,et al.A versatile wavelet domain noisefiltration technique formedical imaging[J].IEEE Transactions on Medical Imaging,2003,22(3):323-331.

[5] PUVANATHASAN P,BIZHEVA K.Interval type-II fuzzy anisotropic diffusion algorithm for speckle noise reduction in optical coherence tomography images[J].Optics Express,2009,17(2):733-746.

[6] 柯麗,杜強(qiáng),蘇哲.應(yīng)用多級(jí)維納濾波的OCT圖像除噪方法[J].光學(xué)精密工程,2008,16(4):740-745.

[7] 蘇鵬,孫延奎,田小林.采用雙樹(shù)復(fù)小波和混合概率模型的光學(xué)相干層析圖像去噪[J].應(yīng)用科學(xué)學(xué)報(bào),2011,29(5):467-472.

[8] MACIEJ S,IWONA G,DANIEL S.Efficient reduction of speckle noise in Optical Coherence Tomography[J].Optics Express,2012,20(2):1337-1359.

[9] SCHMITT J,XIANG S,YUNG K.Speckle in optical coherence tomography[J].Journal of Biomedical Optics,1999,4(1),95-105.

[10] 劉新文,王惠楠,錢志余.小波變換對(duì)OCT圖像的降噪處理[J].光子學(xué)報(bào),2006,35(6):935-939.

[11] 胡鵬.基于光學(xué)相干層析術(shù)的圖像處理研究[D].南京:南京理工大學(xué)碩士學(xué)位論文,2010.

[12] 李自勤,王騏,李琦,等.激光成像雷達(dá)系統(tǒng)中散斑像的乘法模型及其濾除[J].中國(guó)激光,2003,30(8):717-720.

[13] WONG A,MISHRA A,BIZHEVA K,et al.General Bayesian estimation for speckle noise reduction in optical coherence tomography retinal imagery[J].Optics Express,2010,18(8):8338-8352.

[14] 王博.高噪聲率紅外圖像直方圖加權(quán)濾波算法[J].紅外與毫米波學(xué)報(bào),2007,26(5):380-385.

[15] 甘俊英,陳銀河,高建虎.基于加權(quán)直方圖和均值漂移的實(shí)時(shí)跟蹤算法[J].計(jì)算機(jī)仿真,2008,25(11):208-211.

[16] 王曉軍,孫洪.SAR圖像相干斑濾波性能評(píng)估模型[J].雷達(dá)科學(xué)與技術(shù),2008,6(1):33-34.