光學(xué)相干層析成像的視網(wǎng)膜層狀結(jié)構(gòu)自動分割

高用賀，李躍杰,，王立偉，張明蓉

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

2 天津市眼科醫(yī)學(xué)設(shè)備技術(shù)工程中心，天津市，300384

光學(xué)相干層析成像的視網(wǎng)膜層狀結(jié)構(gòu)自動分割

【作 者】高用賀1，李躍杰1,2，王立偉1，張明蓉1

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

2 天津市眼科醫(yī)學(xué)設(shè)備技術(shù)工程中心，天津市，300384

目的 利用算法實(shí)現(xiàn)對視網(wǎng)膜層狀結(jié)構(gòu)的自動分割及定量分析是光學(xué)相干層析成像技術(shù)應(yīng)用于青光眼及視網(wǎng)膜病變早期診斷的關(guān)鍵?，F(xiàn)存處理方法對圖像質(zhì)量要求較高且可靠性不高。該文擬利用改進(jìn)的非線性復(fù)合擴(kuò)散濾波等方法解決這個問題。方法 首先對自主搭建的OCT系統(tǒng)獲得的20幅視網(wǎng)膜圖像，通過自動閾值、改進(jìn)的非線性復(fù)合擴(kuò)散濾波、形態(tài)學(xué)操作、峰值探測等綜合算法，進(jìn)行分割，比較準(zhǔn)確的分割出內(nèi)界膜（ILM）、外核層（ONL）、內(nèi)節(jié)層和外節(jié)層（IS/OS）以及視網(wǎng)膜色素上皮與脈絡(luò)膜層（RPE_ChCap）邊界，最后測量得到視網(wǎng)膜的厚度。結(jié)果 本算法對視網(wǎng)膜的分割與專家手動測量有較好的一致性，視網(wǎng)膜中心凹測量結(jié)果與Zeiss Stratus OCT視網(wǎng)膜中心厚度212±20 μm數(shù)據(jù)一致。結(jié)論該文提出的算法有希望應(yīng)用于臨床視網(wǎng)膜疾病的診斷。

視網(wǎng)膜；光學(xué)相干層析；層狀結(jié)構(gòu)；非線性復(fù)擴(kuò)散濾波；形態(tài)學(xué)操作

0　引言

光學(xué)相干層析成像技術(shù)(Optical Coherence Tomography，OCT) 是近年來迅速發(fā)展起來的一種成像技術(shù)，最早由美國麻省理工學(xué)院的 Huang D等在1991年成功應(yīng)用于眼組織[1-2]。研究證實(shí)，黃斑中心凹光感受器層厚度隨青光眼病程進(jìn)展呈現(xiàn)出動態(tài)、非線性變化的趨勢：在青光眼早期明顯增厚，而在青光眼中晚期又逐漸下降，厚度接近正常水平[3]。同時，黃斑區(qū)視網(wǎng)膜厚度也是診斷及監(jiān)測糖尿病黃斑水腫遺傳性黃斑病變等眼科疾病的重要依據(jù)[4]。對這些有重要意義的數(shù)據(jù)的測量則依賴于對光學(xué)相干層析圖像中層狀結(jié)構(gòu)的精確分割。因此，視網(wǎng)膜層狀結(jié)構(gòu)的定量分析成為眼科光學(xué)相干層析成像的一項核心技術(shù)。在此領(lǐng)域國外機(jī)構(gòu)起步早，2001年，Koozekanani D等[5]發(fā)表了視網(wǎng)膜分割的文獻(xiàn)，此后每年都有很多關(guān)于視網(wǎng)膜層狀結(jié)構(gòu)分割的文獻(xiàn)發(fā)表，并在視網(wǎng)膜層狀結(jié)構(gòu)的分割層數(shù)與精確度上逐步提升[6-8]。國內(nèi)OCT成像研究領(lǐng)域進(jìn)展相對較慢，在對OCT視網(wǎng)膜層狀結(jié)構(gòu)的分割領(lǐng)域研究比較少，可檢索到的中文文獻(xiàn)很少，水平相對國外仍有一定差距[9]。這也嚴(yán)重制約了OCT眼科光學(xué)相干層析成像儀國產(chǎn)化的步伐?，F(xiàn)有的視網(wǎng)膜層狀結(jié)構(gòu)分割算法，大多基于多次幀平均(即對眼底某個位置截面掃描多次)后的圖像進(jìn)行分割。圖像噪聲大幅減少，雖然分割結(jié)果比較理想，但在OCT儀器在對人眼掃描過程中，由于很多原因例如人眼總會有不自主抖動，多次掃描經(jīng)常會有偽差，大大影響算法分割的準(zhǔn)確性。本文提出的視網(wǎng)膜層狀結(jié)構(gòu)的自動算法，是基于一幀OCT視網(wǎng)膜A掃描圖像基礎(chǔ)上，先用改進(jìn)的自適應(yīng)非線性復(fù)合擴(kuò)散濾波器[10]減少噪聲影響，再使用形態(tài)學(xué)操作等方法進(jìn)行分割。用自主研發(fā)的OCT平臺獲得的視網(wǎng)膜圖像進(jìn)行算法測試，對視網(wǎng)膜厚度的測量及層狀結(jié)構(gòu)的定量分析取得了良好的結(jié)果。

1　視網(wǎng)膜層狀結(jié)構(gòu)分割

利用天津市眼科醫(yī)學(xué)設(shè)備技術(shù)工程中心所搭建的高分辨率SD-OCT樣機(jī)進(jìn)行眼底視網(wǎng)膜掃描。醫(yī)學(xué)上，視網(wǎng)膜結(jié)構(gòu)可分為七層，分別為：視網(wǎng)膜神經(jīng)纖維層(retinal never fiber layer, RNFL)、神經(jīng)節(jié)細(xì)胞層和內(nèi)網(wǎng)狀層(inner plexiform layer, IPL)、內(nèi)核層(inner nuclear layer, INL)、外網(wǎng)狀層(outer plexiform layer, OPL)、光感受器內(nèi)節(jié)/外節(jié)(inner and outer photoreceptor segment, IS/OS)、視網(wǎng)膜色素上皮與脈絡(luò)膜層(retinal pigment epithelium and choroid, RPE_ ChCap)，如圖1所示。

圖1　OCT視網(wǎng)膜圖分層結(jié)構(gòu)Fig.1　Layers structures of retina picture.

1.1內(nèi)界膜(Internal limiting membrane, ILM)定位

首先，對原OCT視網(wǎng)膜圖像使用最大類間方差法(簡稱OTSU)找到一個合適的閾值(threshold)，這個閾值在[0, 1]范圍內(nèi)。再利用這個閾值，用Matlab自帶的im2bw函數(shù)將灰度圖像轉(zhuǎn)換為二值圖像。此方法比人為手動設(shè)定的閾值效果更好且方便。對二值圖像水平方向上(即0o)實(shí)施線形膨脹運(yùn)算，結(jié)構(gòu)元素參數(shù)設(shè)為4，這樣可以將ILM邊界填充的更平滑。遍歷尋找圖像每一列值為1的縱坐標(biāo)最小的位置，即ILM初始邊界。對初始邊界第6個點(diǎn)開始，每一個點(diǎn)與其前5個點(diǎn)坐標(biāo)均值作比較，若差的絕對值大于5，則用此坐標(biāo)均值代替此點(diǎn)的坐標(biāo)，以達(dá)到平滑邊界的目的，最終處理結(jié)果為ILM邊界，如圖2所示。

圖2　ILM（內(nèi)界膜）層邊界Fig.2　Inner Limiting Membrane

1.2RPE_Ch和IS/OS邊界定位

1.2.1圖像預(yù)處理

由于OCT系統(tǒng)中，目標(biāo)的光學(xué)性能和移動、光源的光斑的大小和時間相干性、傳輸光的多次散射和相位畸變等都是造成OCT圖像散斑噪聲的主要因素。為了對OCT視網(wǎng)膜圖像更準(zhǔn)確的層分割，顯然在保持圖像層狀信息的基礎(chǔ)上，去除影響圖像質(zhì)量的散斑噪聲，顯得尤為重要。

目前，已有一些算法，如中值濾波，高斯濾波，小波非線性閾值濾波，零幅度算法(Zero-adjustment procedure, ZAP)[11]，清潔算法(CLEAN)[12-13]，增強(qiáng)的Lee濾波器(ELEE)等方法。本文采用一種改進(jìn)的自適應(yīng)非線性復(fù)擴(kuò)散濾波器，是由Rui Bernardes于2010年提出的，它基于普通的非線性復(fù)合擴(kuò)散濾波器。相對于其他算法，此算法對OCT視網(wǎng)膜圖像具有更好的噪聲消除和層狀結(jié)構(gòu)邊緣保持能力。

一般的非線性各向異性復(fù)合擴(kuò)散是基于式(1)的解：

式中，△●是散度，△是梯度，D是擴(kuò)散系數(shù)。

將式（1）時域前向離散并化為空域有限差分式，即

式中，△和△分別為離散拉普拉斯和梯度算

hh子，是迭代步數(shù)n的時間，i, j, m是體數(shù)據(jù)塊的坐標(biāo)向量，在2D數(shù)據(jù)中

成立，在（4）式中我們令△h=△y=1且對于2D圖像數(shù)據(jù)α=4。Salinas將擴(kuò)散系數(shù)定義為

式中，i =√-1，k是閾值參數(shù)，θ是一個接近于0的相位角。此算法中，擴(kuò)散使得圖像平滑但衰減了層狀結(jié)構(gòu)邊緣信息。

擴(kuò)散系數(shù)可以簡化為

基于OCT視網(wǎng)膜圖像的背景（主要是玻璃體、房水）區(qū)域強(qiáng)度較低，而視網(wǎng)膜區(qū)域強(qiáng)度相對較高的基礎(chǔ)上，在低強(qiáng)度區(qū)域促進(jìn)擴(kuò)散，而在高強(qiáng)度區(qū)域降低擴(kuò)散，這樣能更好的保留層狀結(jié)構(gòu)邊緣信息。因此，要改進(jìn)參數(shù)k，使其在圖像低強(qiáng)度區(qū)域增強(qiáng)，在圖像高強(qiáng)度區(qū)域減小，如下式

式中，g=GN,σ*Re(I)，GN,σ是高斯平均內(nèi)核，*是卷積運(yùn)算，尺寸為N×N或N×N×N，σ為標(biāo)準(zhǔn)差。

相對于其他非線性復(fù)合擴(kuò)散算法，大多設(shè)置的是固定的時間步長(△t)，由于擴(kuò)散系數(shù)依賴于圖像的梯度，并且在擴(kuò)散迭代初始時因?yàn)樵肼暤木壒适固荻戎蹈?，故這里設(shè)置為自適應(yīng)時間步長，即

式中，α=4，且a+b≤1。

將擴(kuò)散時間最大值設(shè)為0.8 s，對OCT視網(wǎng)膜圖像進(jìn)行處理，之后截取ILM層上部一塊兒矩形背景區(qū)域，求取平均像素強(qiáng)度Ibg，用下式對OCT視網(wǎng)膜圖像進(jìn)行局部增強(qiáng)，即

式中，I為整體OCT視網(wǎng)膜圖像，Ibg是所截取部分的平均強(qiáng)度，這里β=2，重復(fù)連續(xù)進(jìn)行3次。以達(dá)到增強(qiáng)RPE部分的目的，結(jié)果如圖4所示。

1.2.2形態(tài)學(xué)分割RPE_ChCap

二值形態(tài)學(xué)中的運(yùn)算對象是集合，通常給出一個圖像集合和一個結(jié)構(gòu)元素集合，利用結(jié)構(gòu)元素對圖像進(jìn)行操作。其中，結(jié)構(gòu)元素是一個用來定義形態(tài)操作中所用到的形狀和大小的矩陣，該矩陣僅由0和1組成，可以具有任意的大小和維數(shù)，數(shù)值1代表鄰域內(nèi)的像數(shù)，形態(tài)學(xué)運(yùn)算都是對數(shù)值為1的區(qū)域進(jìn)行運(yùn)算。

圖4　增強(qiáng)RPE結(jié)果Fig.4　Results of enhanced RPE

形態(tài)學(xué)開運(yùn)算，即對其先腐蝕后膨脹?？上?xì)小物體、在纖細(xì)處分離物體、平滑較大物體的邊界時不明顯改變其面積。

首先，同ILM層邊界定位中一樣，將圖4二值化，用Matlab中bwareaopen函數(shù)去掉面積70以下區(qū)域，結(jié)果如圖5所示。

圖5　二值化結(jié)果Fig.5　Results of binarization

接下來開運(yùn)算操作，為了使OCT視網(wǎng)膜層狀結(jié)構(gòu)邊緣更加準(zhǔn)確，腐蝕操作結(jié)構(gòu)元素選用半徑為2的圓盤狀結(jié)構(gòu)；膨脹操作結(jié)構(gòu)元素選用長度為20，方向?yàn)樗剑?o）的線形結(jié)構(gòu)，用Matlab中bwareaopen函數(shù)去掉面積5 000以下區(qū)域（即RNFL區(qū)域）。遍歷尋找圖像每一列值為1的縱坐標(biāo)最小、最大的位置，即分別為IS/OS邊界、RPE_Ch邊界。為了減小IS/OS層邊界左右兩端由于線型膨脹操作所引起的誤差，以RPE_Ch邊界最低點(diǎn)為基準(zhǔn)，將圖像逐列下移，使RPE_Ch邊界拉平為水平直線。

對拉平的圖像以相同設(shè)置重復(fù)如上操作，即對RPE部分局部增強(qiáng)、形態(tài)學(xué)開運(yùn)算分割，并邊界定位，結(jié)果如圖6所示。

圖6　拉平原圖后RPE層分割結(jié)果Fig.6　Results of segmentation of RPE after alignment

1.2.3分割ONL

對經(jīng)過預(yù)處理過的拉平的OCT視網(wǎng)膜圖像每11列像素寬的A掃描疊加求均值賦于第1列，以第200列像素為例，畫出其灰度強(qiáng)度圖。最大的峰值出現(xiàn)在RPE_ChCap層內(nèi)部，從最大峰值處向上使用峰值探測法尋找，第一個峰值即是ONL層內(nèi)邊界。之后，遍歷圖像所有列，用上述方法找到所有峰值點(diǎn)。將所有峰值點(diǎn)連成光滑曲線，即為ONL層邊界，如圖7所示，由上往下第二條分割線即為ONL層邊界。

圖7　ONL層位置Fig.7　Location of ONL

2　視網(wǎng)膜厚度測量

將ILM邊界距離RPE_ChCap邊界的兩個極大值分別定義為旁中心凹（parafovea），在旁中心凹中間尋找ILM邊界的最低點(diǎn)并定義為中心凹（fovea）如圖8所示。提取ILM邊界與IS/OS邊界的像素差，再經(jīng)過定標(biāo)后生成視網(wǎng)膜層厚度分布圖，如圖9所示。對所有列厚度加權(quán)求平均值，即為此A掃描的視網(wǎng)膜平均厚度。視網(wǎng)膜光感受器層厚度是以中心凹處ILM邊界到IS/OS邊界處的距離表示，如圖8箭頭長度所示。

圖8　中心凹及旁中心凹位置Fig.8　Location of fovea and parafove

本文在檢驗(yàn)算法的準(zhǔn)確性及可靠性的過程中，用自主搭建的設(shè)備采集20個正常人眼底視網(wǎng)膜OCT圖像。運(yùn)用本文所提出的算法進(jìn)行自動測量黃斑中心凹光感受器厚度和視網(wǎng)膜平均厚度。并與專家手動分割測量結(jié)果進(jìn)行比較，結(jié)果如表1所示。

由表1可知，本文提出的自動算法與專家手動分割測量結(jié)果無顯著性差異，具有較好的一致性。其中，視網(wǎng)膜中心凹厚度與Zeiss Stratus OCT視網(wǎng)膜中心厚度212±20 μm數(shù)據(jù)一致[14]。之所以自動測量的結(jié)果較手動測量偏低，是由于形態(tài)學(xué)操作的誤差所致，但差異范圍在3%以內(nèi)。

圖9　視網(wǎng)膜厚度分布曲線Fig.9　Retinal thickness distribution curve

表1　黃斑中心凹光感受器厚度和視網(wǎng)膜平均厚度(n=20) (μm)Tab.1　Photoreceptor thickness of macular fovea and mean retinal thickness(n=20)

3　結(jié)論

光學(xué)相干層析系統(tǒng)具有非接觸、無損傷、實(shí)時、高分辨力等特性，對光學(xué)相干層析圖像中的層狀結(jié)構(gòu)進(jìn)行定量提取與其他檢測方法相比具有很多優(yōu)點(diǎn)。本文提出的算法，首先通過使用最大類間方差法設(shè)定閾值結(jié)合線型膨脹運(yùn)算對ILM層邊界定位。其次，通過改進(jìn)的自適應(yīng)非線性復(fù)合擴(kuò)散濾波器對OCT視網(wǎng)膜圖像進(jìn)行預(yù)處理，既濾除圖像中的噪聲，同時又保留了更多的圖像層狀特征，減少了層狀結(jié)構(gòu)位置失真。對OCT視網(wǎng)膜圖像的預(yù)處理，使得背景區(qū)域灰度更加均勻集中，對RPE_ChCap區(qū)域進(jìn)行線性增強(qiáng)奠定了基礎(chǔ)，進(jìn)而對IS/OS邊界和RPE_ChCap邊界進(jìn)行定位。最后，利用峰值探測法對ONL層進(jìn)行邊界定位。結(jié)果顯示，此方法具有較好的自適應(yīng)性和重復(fù)性，可以應(yīng)用在臨床OCT系統(tǒng)中。

當(dāng)然，此方法也有需要改進(jìn)的地方，還不能對視網(wǎng)膜全部層狀結(jié)構(gòu)進(jìn)行分割。對較嚴(yán)重眼底病變的視網(wǎng)膜圖像自動分層的準(zhǔn)確度有待提高。今后，可在此算法基礎(chǔ)上，利用馬爾可夫隨機(jī)場模型對視網(wǎng)膜OCT圖像的全部層狀結(jié)構(gòu)進(jìn)行分割，并添加模式識別，以實(shí)現(xiàn)視網(wǎng)膜病變的全自動識別診斷。

[1] Huang D, Swanson EA, Lin CP, et al. Optical coherence tomography[J]. Science, 1991, 254(5035): 1178-1181.

[2] Drexler W, Fujimoto JG. State-of-the-art retinal optical coherence tomography[J]. Prog Retin Eye Res, 2008, 27(1): 45-48.

[3] Quigley HA, Nickells RW, Kerrigan LA, et a1．Retinal ganglion cell death in experimental glaucoma and after axotomy occurs by apoptosis[J]. Invest Ophthalmol Vis Sci, 1995, 36: 774-786．

[4] Ishikawa H, Stein DM, Wollstein G, et al. Macular Segmentation with Optical Coherence Tomography[J]. Invest Ophthalmol Vis Sci, 2005, 46(6): 2012-2017.

[5] Koozekanani D, Boyer K, Roberts C. Retinal thickness measurements from optical coherence tomography using a Markov boundary model[J]. IEEE Trans Med Imag, 2001, 20(9): 900-916.

[6] Tan O, Chopra V, Lu ATH, et al. Detection of macular ganglion cell loss in glaucoma by Fourier-domain optical coherence tomography[J]. Ophthalmology, 2009, 116(12), 2305-2314.

[7] Yazdanpanah A, Hamarneh G, Smith B, et al. Intra-retinal layer segmentation in optical coherence tomography using an active contour approach[C]. Proc MICCAI, 2009: 649-656.

[8] Kajic V, Povazay B, Hermann B, et al. Robust segmentation of intraretinal layers in the normal human fovea using a novel statistical model based on texture and shape analysis[J]. Opt Express, 2010, 18(14): 14730-14744.

[9] 江源源, 周傳清, 任秋實(shí). 基于光學(xué)相干層析的視網(wǎng)膜圖像分割[J]. 北京生物醫(yī)學(xué)工程, 2011, 30(5): 452-456.

[10] Bernardes R, Maduro C, Serranho P, et al. Improved adaptive complex diffusion despeckling fi lter[J]. Opt Express, 2010, 18(23): 24048-24059.

[11] Yung KM, Lee SL, Schmitt JM. Phase-domain processing of optical coherence tomography images[J]. J Biomed Opt, 1999, 4(1): 125-136.

[12] Fried DL. Analysis of the clean algorithm and implications for super resolution[J]. J Opt Soc Am A,1995, 12(5): 853-860.

[13] Schmitt JM. Restoration of optical coherence images of living tissue using the clean algorithm[J]. J Biomed Opt, 1998, 3(1): 66-75.

[14] Legarreta JE, Gregori G, Punjani OS, et al. Macular thickness measurements in normal eyes using spectral domain optical coherence tomography[J]. Ophthalmic Surg Las Imag, 2008, 39(4 Suppl): S43-49.

Automated Segmentation of Retina Layer Structures on Optical Coherence Tomography

【 Writers 】Gao Yonghe1, Li Yuejie1,2, Wang Liwei1, Zhang Mingrong1
1 Peking Union Medical College&Chinese Academy of Medical Sciences Institute of Biomedical Engineering, Tianjin, 300192
2 Tianjin Ophthalmic Medical Device Technology Engineering Center, Tianjin, 300384

retina, optical coherence tomography, layer structure, nonlinear complex diffusion fi ltering, morphological operations

TP391.41

A

10.3969/j.issn.1671-7104.2014.02.004

1671-7104(2014)02-0094-04

2014-01-02

中國醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)工程研究所院所基金（1312）；天津科技創(chuàng)新體系及平臺建設(shè)計劃項目（13ZXGCCX06300）

高用賀，研究生，E-mail: gyhklgyhkl@126.com

李躍杰，研究員，E-mail: liyj_1@sina.com

【 Abstract 】Objective Using the algorithm on the layered structure of the retina and quantitative analysis of the automatic segmentation technique is the key to the early diagnosis of glaucoma and other retinopathy on optical coherence tomography. Existing methods require high qulity image and have low reliability. This paper used the improved complex nonlinear diffuse fi ltering and other methods to solve this problem. Methods This paper includes algorithm such as automatic threshold, improved complex nonlinear diffusion filtering, morphological operations and peak detection. Use the method for the segmentation of 20 retinal layers images which acquired on the self-builded OCT system, the boundary of inner limiting membrane(ILM), outer nuclear layer(ONL), the photoreceptor segments(IS/ OS) and the RPE_ChCap layer are detected accurately. At last, the photoreceptor layer thickness is measured. Results The results of segmentation and measurement are good corresponded with expert manual segmentation and measurements, retinal foveal measurements data is consistent with Zeiss Stratus OCT central retinal thickness 212±20 μm. Conclusion The algorithm proposed is prospective applied to clinical diagnosis of retinal diseases.