多算法融合的視盤分割方法

傅迎華，王雅靜，付東翔，楊振宇

(上海理工大學(xué) 光電信息與計算機工程學(xué)院， 上海 200093)

1 引 言

在正常視網(wǎng)膜圖像中，視盤一般表現(xiàn)為一個近似圓形、亮黃色的物體.對視網(wǎng)膜眼底圖像中的視盤區(qū)域進(jìn)行定位和分割，一方面能夠為眼底疾病如青光眼的診斷提供必要的輔助信息；另一方面，由于一些眼底病變同視盤具有類似的亮度信息，以及排除某些病變的發(fā)生區(qū)域，對其進(jìn)行定位和分割，能夠有效地降低這些相關(guān)病變的特異性.

對于視盤分割的方法根據(jù)視盤的結(jié)構(gòu)特征，大致分為4類：

1)基于邊緣的方法：根據(jù)視盤的灰度信息來提取視盤邊界，獲得視盤輪廓達(dá)到分割的目的.Xu等[1]提出了一種可變形模型，將主動輪廓形變與聚類相結(jié)合.Lowell等[2]利用全局橢圓參數(shù)模型結(jié)合局部變形模型分割OD邊界.該方法對初始曲線位置較為敏感，且在低對比度圖像中收斂效果不明顯.

2)基于形狀的方法：根據(jù)視盤類圓的形狀先驗利用不同的算法達(dá)到視盤分割的目的.Morales等[3]使用橢圓或圓來擬合OD邊界.該類方法的優(yōu)點是簡單直觀，但其前提是默認(rèn)視盤為圓形.然而視盤往往不是一個規(guī)則的圓形或橢圓形，因此往往會漏掉實際的視盤邊界.

3)基于區(qū)域的方法：將眼底圖像分割成不同的區(qū)域，再根據(jù)視盤的特點找出視盤所在的區(qū)域達(dá)到分割的目的.Reza等提出的基于控制標(biāo)記的分水嶺變換算法[4]，通過最小值變換的方法得到內(nèi)部標(biāo)記并對梯度圖像進(jìn)行修改以控制過分割.該類算法易受明亮病灶和血管邊緣反光的干擾，魯棒性并不高.

4)由于以上三類方法各有優(yōu)缺點，選取適當(dāng)?shù)姆椒ㄟM(jìn)行算法組合有時會獲得更好的效果.Abdullah[5]等提出將形態(tài)學(xué)運算、霍夫圓變換和生長切割算法相結(jié)合分別用于視盤的增強、定位和分割，該方法能夠準(zhǔn)確地定位視盤，視盤的分割效果在部分?jǐn)?shù)據(jù)集上有明顯優(yōu)勢.曾杰[6]提出用基于梯度差與DOG算子的方法定位視盤，提取視盤的局部候選區(qū)域再利用基于標(biāo)記的分水嶺變換進(jìn)行區(qū)域篩選，最終利用圓擬合和活動輪廓模型提取視盤邊緣.Qureshi[7]等提出融合金字塔分解、邊緣檢測、熵濾波、霍夫變換和基于特征向量和均勻采樣網(wǎng)格5種算法建立置信度來檢測視盤.該類方法將不同原理的算法組合利用算法的優(yōu)勢彌補單個算法的不足，具有良好的檢測效果.

圖1 視盤分割的流程Fig.1 Flow chart of optic disc segmentation

本文融合多種傳統(tǒng)的分割方法對視網(wǎng)膜眼底圖像進(jìn)行視盤分割：對于分辨率為m×n的彩色眼底圖像，預(yù)處理部分根據(jù)主成分分析提取主要的視盤信息作為其灰度圖片，利用形態(tài)學(xué)開閉操作去除血管.接著先用霍夫圓變換的方法對視盤進(jìn)行粗定位得到視盤中心，再用分水嶺變換分割視盤，并將得到的視盤輪廓作為初始曲線.最后利用CV模型將初始曲線演化到真正的視盤邊界得到最終的分割結(jié)果，算法流程如圖1所示.由于分水嶺算法會產(chǎn)生過分割導(dǎo)致分割結(jié)果不正確，本文采用改進(jìn)的基于霍夫圓變換的方法來彌補過分割的不足得到準(zhǔn)確的視盤邊界.

2 預(yù)處理

由于視盤是視網(wǎng)膜主血管的匯集區(qū)域，視盤的分割受血管的影響較大，很難進(jìn)行準(zhǔn)確定位與分割.因此，需要在保留視盤信息的前提下去除血管.

2.1 基于PCA的圖像灰度化

考慮到視盤的主要信息體現(xiàn)在紅色通道(亮度特征)和綠色通道(信息完整)中，本文采用主成分分析(Principal Components Analysis，PCA)的方法對圖像進(jìn)行灰度化處理.首先，將R,G,B三個顏色空間的數(shù)據(jù)作均值為0的歸一化處理；再求其協(xié)方差矩陣，如公式(1)所示；接著計算協(xié)方差矩陣C的特征值λ和特征向量u；最后將三個顏色空間的值投影到特征值最大的特征向量上，完成圖像的灰度化[6].

(1)

PCA的方法綜合了三個通道有意義的信息，第一成分保留了完整的視盤結(jié)構(gòu)信息，且血管信息被抑制；第二成分血管信息較為明顯，但對比度高，突顯了視盤的明亮信息；第三成分中視盤信息缺失嚴(yán)重，明亮病變保留完整.如圖2中(a)所示，為了保證圖像保留絕大部分的結(jié)構(gòu)和對比度特征，本文選取視盤信息較明顯的第一成分I1為主成分結(jié)合第二成分I2作為灰度圖片，如公式(2)所示.

IPCA=0.8*I1+0.2*I2

(2)

2.2 形態(tài)學(xué)去除血管

本文采用基于形態(tài)學(xué)的開運算和閉運算方法去除血管，以減少血管對于視盤分割的影響，如公式(3)和公式(4)所示.采用先閉后開的順序消除血管并平滑圖像，如公式(5)所示.

形態(tài)學(xué)開運算：I0=IPCA°k=(IPCA?k)⊕k

(3)

形態(tài)學(xué)閉運算：IC=IPCA·k=(IPCA⊕k)?k

(4)

先閉后開：IRV=(IC?k)⊕k

(5)

其中，IPCA表示PCA處理后得到的灰度圖像，k表示用于相關(guān)操作的核，IRV表示血管消除后的圖像.

圖2 消除血管結(jié)果對比Fig.2 Comparison of the vessel removed result

圖2(b)表明基于先閉后開的形態(tài)學(xué)方法能夠有效地去除掉血管及圖像中小的病灶.

3 霍夫圓變換定位視盤

由于視盤在眼底圖像中表現(xiàn)為明亮的類圓，而圓心半徑均未知，本文采用改進(jìn)的基于霍夫變換的圓檢測方法能夠有效檢測出圖像中未知圓心和半徑的圓形物體.

3.1 霍夫圓變換

在圓心半徑均未知的情況下需要預(yù)先估計檢測目標(biāo)(視盤)半徑的范圍(rmin,rmax)，再根據(jù)圖像平面上的點進(jìn)行求解[8].先將灰度圖像做閾值分割，如公式(6)，再進(jìn)行邊緣檢測提取視盤骨架作為檢測對象.這一步可以有效地去除掉其他干擾信息.然后，根據(jù)圓的一般方程變形形式進(jìn)行Hough變換，將圖像空間轉(zhuǎn)換成參數(shù)空間對所有可能落在邊緣上的點進(jìn)行統(tǒng)計計算找出視盤的中心和半徑，如公式(7).

(6)

(7)

其中，T是閾值，IS為閾值分割后的圖像，(a,b)表示圓心，r是半徑，(x,y)是圖像平面上的點.對于視盤的圖像元素IS(i,j)=1，對于背景的圖像元素IS(i,j)=0.

3.2 改進(jìn)的霍夫圓變換

已經(jīng)去除血管的圖像IRV中的視盤輪廓較明顯且灰度值較大，也就是說視盤區(qū)域較亮而背景較暗，根據(jù)這一特性本文通過直方圖抽取灰度級的方法自動地尋找合適的閾值，提取出準(zhǔn)確的視盤骨架從而提高視盤定位的精確度[9]. 如圖3所示，根據(jù)直方圖公式將視盤區(qū)域所在的較大灰度級求均值得到閾值T，如公式(8)所示.

圖3 直方圖選取閾值Fig.3 Histogram selection threshold

(8)

其中，rj表示第j級灰度值，nj表示圖像中灰度為rj的像素個數(shù)，L是圖像中可能的灰度級數(shù)量，G表示方差較小的較大灰度級.

基于霍夫變換的圓檢測方法能夠有效檢測出圖像中未知圓心和半徑的圓形物體.但也正由于圓心和半徑均未知，通常一個目標(biāo)檢測圖像會檢測出不止一個圓.本文在霍夫圓變換檢測圓的基礎(chǔ)上，根據(jù)n個檢測圓分別與檢測對象相交區(qū)域面積的大小作為評價指標(biāo)來確定真實的視盤輪廓及圓心，如公式(9).

(9)

圖4 視盤定位Fig.4 Optic disc location

其中，D表示霍夫圓變換檢測出的圓,O是檢測對象,Area()是求面積的函數(shù).Si的值越大表明檢測出的圓與檢測對象的重疊率最高.Simax對應(yīng)的檢測圓的圓心就是視盤中心，如圖4所示. (a)為Hough圓檢測的視盤定位結(jié)果，(b)中的“+”號為視盤中心，虛線為根據(jù)檢測出的視盤半徑所作的圓.

4 改進(jìn)的傳統(tǒng)的視盤分割算法

4.1 改進(jìn)的分水嶺分割算法

基于控制標(biāo)記的分水嶺算法大多利用局部極小值的方法獲取內(nèi)部標(biāo)記[10].算法從參考圖像的全局極小值點開始，小于或等于某一平面的集水盆地(如公式(10)所示)會分配一個唯一的標(biāo)記，而高于這一平面的元素，若其鄰域內(nèi)已有標(biāo)記，則給它分配相同的標(biāo)記；若其鄰域內(nèi)無標(biāo)記，則給它分配一個新的標(biāo)記.進(jìn)行反復(fù)標(biāo)記，直到圖像中的每個像素都處于某個集水盆地為止，分水嶺變換則完成[11].對于視網(wǎng)膜眼底圖像，其視盤內(nèi)部區(qū)域?qū)?yīng)集水盆地的最低點，視盤邊緣對應(yīng)分水嶺.但是由于視網(wǎng)膜眼底圖像中明亮病變和血管邊緣的反光的存在，導(dǎo)致一幅眼底圖像中存在多個局部極值，直接使用分水嶺算法很容易產(chǎn)生過分割現(xiàn)象，如圖5 (a)所示.

T[m]={(s,t)|GRV(s,t)

(10)

幾何意義上，T[m]是梯度圖像GRV(x,y)中位于平面GRV(x,y)=m下方的點的坐標(biāo)集合.梯度閾值m從圖像梯度的最小值開始依次取整遞增，min和max分別表示圖像梯度范圍的最小值和最大值.

為了避免過分割現(xiàn)象的產(chǎn)生，本文提出了基于霍夫圓變換的分水嶺算法，即在霍夫圓變換定位視盤的基礎(chǔ)上，將視盤的中心坐標(biāo)作為標(biāo)記lin，并應(yīng)用為局部最低值.將視盤定位得到的圓形輪廓擴大1.5倍作為外部標(biāo)記lout，得到最終的控制標(biāo)記如公式(11)所示.

l=lin∪lout

(11)

圖5 視盤分割結(jié)果Fig.5 Optic disc segmentation results

圖5 (b)為本文算法的分割結(jié)果，圖5 (c)-圖5(d)中粗實線表明當(dāng)眼底圖像中視杯過于明亮?xí)r，基于霍夫變換的分水嶺分割算法可能會導(dǎo)致糟糕的分割結(jié)果，或者有些視盤的邊緣受血管的影響并沒有被檢測出來.為了避免這一問題，接下來本文采用了改進(jìn)的活動輪廓模型對該算法得到的曲線輪廓進(jìn)行修正，將其逼近到視盤真正的邊緣，提高算法的魯棒性.

4.2 CV模型

Chan-Vese模型[12]是基于區(qū)域的水平集分割算法，以待分割圖像的目標(biāo)區(qū)域和背景區(qū)域的像素灰度平均值差異最小為準(zhǔn)則來分割圖像，其能量泛函定義為：

(12)

其中，C表示演化曲線，μ為曲線C的長度權(quán)值，v為曲線C包含區(qū)域的面積權(quán)值，μ,v≥0.λ1，λ1分別為曲線C內(nèi)外部區(qū)域能量函數(shù)的權(quán)值，λ1，λ1>0.in(c)和out(c)分別對應(yīng)曲線C的內(nèi)部和外部區(qū)域.c1和c2分別表示圖像中分割物體內(nèi)部和外部的亮度均值.當(dāng)曲線C恰好為物體邊界時，能量ECV(c1,c2,C)最小.

傳統(tǒng)的CV模型存在收斂速度較慢的問題，并且在目標(biāo)區(qū)域與背景區(qū)域?qū)Ρ榷鹊偷那闆r下算法收斂效果不好.文獻(xiàn)13中提出了一種改進(jìn)的CV模型，該算法加入了自適應(yīng)權(quán)值δ1和δ2，提高目標(biāo)區(qū)域和背景區(qū)域擬合中心的準(zhǔn)確性來改善分割效果，并對λ1，λ2進(jìn)行優(yōu)化加速模型演化[13].改進(jìn)的能量泛函對應(yīng)Heaviside函數(shù)的變形為：

ECV(b1，b2,β1,β2,C)=μ?ΩH′(φ)|▽φ|dxdy+v?ΩH(φ)dxdy+β1?Ω|IRV(x,y)-b1|2H(φ)dxdy+β2?Ω|IRV(x,y)-b2|2(1-H(φ)dxdy

(13)

1http://www.ia.uned.es/～ejcarmona/DRIONS-DB.html.

其中，φ表示0水平集曲線，Ω表示圖像的整體區(qū)域，H(φ)表示理想階躍函數(shù)，如公式(14)所示.H′(φ)表示其一階微分，如公式(15)所示.b1，b2分別為引入自適應(yīng)δ1和δ2后的擬合中心，如公式(16),公式(17)所示.β1，β2表示優(yōu)化的曲線C內(nèi)外區(qū)域能量函數(shù)權(quán)值，為區(qū)域像素灰度的中值絕對差.

(14)

(15)

式中，ε表示趨近于0的正數(shù).

(16)

式中，δ1和δ2定義為：

(17)

式中，k為常數(shù)，避免出現(xiàn)分母為零的情況.

圖5 (c)-圖5(d)中虛線為用改進(jìn)的CV模型將分水嶺變換得到的視盤輪廓推動到真正的視盤輪廓的結(jié)果，與ground truth(細(xì)實線)對比發(fā)現(xiàn)該算法彌補了受視杯過于明亮和血管的影響產(chǎn)生的不足，提高了算法的魯棒性.

5 實驗結(jié)果與分析

本文選取兩個公共數(shù)據(jù)庫對本文方法進(jìn)行實驗和評價，分別為：MESSIDOR[14]和DRIONS-DB1.MESSIDOR數(shù)據(jù)庫共有1200幅具有不同分別率的眼底圖像，包含普通圖像和病灶圖像.為方便實驗，文本使用的眼底圖像均調(diào)整為分辨率為500×500 像素的去光照圖像，所有圖像的視盤邊界由專家手動劃定.DRIONS-DB數(shù)據(jù)庫共有110張分辨率為600×400像素的眼底圖像，其中部分圖像存在光偽影、視乳頭萎縮、邊緣模糊等干擾.該數(shù)據(jù)庫的所有圖像有兩位專家對顯著的視盤輪廓點進(jìn)行標(biāo)記，本文中選取專家1的標(biāo)記結(jié)果作為ground truth.

5.1 視盤定位的實驗結(jié)果

由于視網(wǎng)膜眼底圖像數(shù)據(jù)庫并未提供視盤正確定位的ground truth，本文選取視盤的中心作為ground truth來判斷視盤是否正確定位.通過此評價方法得到本文方法在兩個數(shù)據(jù)庫上的準(zhǔn)確度分別為98%和100%，圖6中“+”號為本文方法的視盤定位結(jié)果.

5.2 視盤分割的實驗結(jié)果

由于MESSIDOR和DRIONS-DB數(shù)據(jù)庫均提供了專家手動標(biāo)記的視盤邊界，因此通過計算重疊率的方法(如公式18)來評估本文方法的視盤分割結(jié)果.

(18)

其中Ap和Ag分別為本文分割的結(jié)果和ground truth，OP的值在0至1之間.OP的值越大，表明分割的結(jié)果與ground truth之間的重疊率越高，即算法的精度越高.如圖6所示，(a)-(b)中虛線為本文方法，細(xì)實線為ground truth；(c)中實線為本文方法，點劃線為ground truth.

圖6 視盤分割結(jié)果Fig.6 Optic disc segmentation results

表1中的數(shù)據(jù)表示在不同重疊率下的圖像所占總圖像的百分比.在DRIONS-DB數(shù)據(jù)庫上的分割結(jié)果明顯優(yōu)于在MESSIDOR數(shù)據(jù)庫，平均重疊率分別為0.94和0.92，因此本文方法具有較好的分割效果.

表1 本文方法在兩個數(shù)據(jù)庫上的重疊率Table 1 Overlap scores of the method in this article on both databases

表2 不同方法在MESSIDOR上的平均重疊率Table 2 Average overlap scores on MESSIDOR for different methods

表2給出了不同方法在MESSIDOR數(shù)據(jù)庫上的平均重疊率.在這些方法中，文獻(xiàn)[3]和文獻(xiàn)[15]使用圓形或橢圓形來擬合視盤邊界[3,15]，由于視盤并非完全規(guī)則的圓形或橢圓形，因此該類方法的重疊率較小，分別為0.82和0.86.文獻(xiàn)[4]通過最小值變換的方法獲得標(biāo)記再進(jìn)行分水嶺變換得到視盤邊界[4]，由于明亮病灶和血管邊緣反光的干擾導(dǎo)致標(biāo)記獲取有誤，該方法的重疊率為0.84，本文改進(jìn)的方法比該方法準(zhǔn)確率提高了8%.文獻(xiàn)[5]和文獻(xiàn)[6]采用了組合算法進(jìn)行視盤檢測和分割，算法的選取和利用各有優(yōu)缺點，其重疊率分別為0.88和0.82.文獻(xiàn)[16]和文獻(xiàn)[17]綜合利用了視盤的邊緣和形狀特點來分割視盤邊界[16,17]，取得了較高的重疊率，分別為0.88和0.89.

6 結(jié) 論

本文提出了融合多種傳統(tǒng)分割算法的視盤分割方法，將霍夫圓變換、分水嶺分割算法和CV模型相結(jié)合.首先通過預(yù)處理和形態(tài)學(xué)的方法去除血管；再用霍夫圓變換的方法定位視盤得到視盤的中心坐標(biāo)；然后基于視盤的中心坐標(biāo)進(jìn)行分水嶺變換得到視盤的初始輪廓；最后利用CV模型將視盤的輪廓進(jìn)行修正，得到最后的視盤邊界.在MESSIDOR和DRIONS-DB數(shù)據(jù)庫上的實驗結(jié)果表明該算法具有很好的視盤定位和分割效果.