基于多尺度濾波的視網膜血管分割算法

梁禮明 盛校棋 藍智敏 錢艷群

(江西理工大學電氣工程與自動化學院 江西 贛州 341000)

0 引 言

血管是視網膜最重要的組成部分之一，視網膜血管分割和血管形態屬性的劃分，如長度、寬度、迂曲和角度，可用于各種心血管和眼科疾病的診斷、篩選、治療和評估[1]，如糖尿病、青光眼[2]、硬性滲出物[3]和高血壓[4]等全身性疾病。然而，通過眼科專家分割視網膜血管存在時效性差，且缺乏準確性。因此，研究血管自動分割和自動識別血管形態屬性的算法成為輔助醫療診斷至關重要的一步。

現有的視網膜分割方法大致可劃分為無監督與有監督方法。其中，無監督學習視網膜血管分割方法可以分為匹配濾波、數學形態學、血管追蹤以及聚類等。文獻[5]提出一種改進的Top-Hat變換應用于視網膜血管分割，使用不同半徑的圓形結構元素來檢測不同寬度的血管，進一步提高微小血管的分割率，但對噪聲較為敏感。文獻[6]提出一種基于貝葉斯概率的血管追蹤方法，結合血管的連通性與灰度級等信息擬合血管結構，該方法較好地解決血管交叉處分割斷裂的問題，但對微小血管分割精度較低。文獻[7]利用主成分分析(PCA)提取血管的特征，并利用閾值進行分割，存在微小血管分割低、血管交叉處分割斷裂等現象。

有監督學習需要專家提供具有標記的金標準視網膜圖像，使用基于局部或全局圖像的一組特征來訓練分類器，充當先驗知識并指導訓練。文獻[8]提出了一種融合多特征的隨機森林分類器，很好地解決了血管交叉處分割斷裂的情況，但存在將視盤誤分割為血管的現象。文獻[9]融合相位特征的血管分割算法，較好地解決了血管相位一致性特征檢測不足的問題，但仍存在微血管分割不足的問題。文獻[10]將血管分割看作為二分類問題，提出一種基于支持向量機(SVM)的混合5D特征區分血管像素與非血管像素，對噪聲有較強的適應性，解決了視盤被誤分割為血管的問題，但存在微小血管易斷裂問題。

因此，針對現有的方法的不足與挑戰，本文提出一種基于多尺度濾波器與AdaBoost分類器相結合的有監督的視網膜分割算法。首先利用二維K-L變換綜合分析圖像的三通道顏色特征信息，解決單一通道存在偽影和噪音大的問題，并采用受限對比度直方圖均衡化(Contrast Limited Adaptive Histogram Equalization, CLAHE)增強視網膜圖像并且移除低頻的噪音，利用雙邊濾波Retinex進行偽影矯正。其次，利用多尺度形態學濾波、多尺度匹配濾波、2D-Gabor濾波及Hessian的各項異性提取較寬和微小血管特征向量，將以上先驗信息與血管特征作為訓練樣本并歸一化。使用專家分割結果作為標簽，通過AdaBoost進行分類，去掉非血管像素，得到初步分割結果。利用AdaBoost分類器與多尺度多類型濾波器相結合，有利于改善分割效果，每個濾波器以不同的方式響應圖像中的不同像素。最后，利用連通域面積等操作去除病灶與偽影，獲得最終的視網膜血管分割結果?？偭鞒倘鐖D1所示。

圖1 本文算法總流程圖

1 數據集與圖像預處理

1.1 數據集

本文采用STARE[11]和DRIVE[12]眼底圖像數據集作為實驗數據。其中：DRIVE數據集各有20張訓練集和測試集；STARE數據集僅含20張眼底圖像，但并沒有劃分訓練集與測試集。故本文利用留一法，選擇一幅圖像作為測試集，并使用數據集中的其他圖像訓練分類器，通過改變測試圖像重復該過程，直到數據集中所有圖像被用于測試一次。

1.2 二維K-L變換

采用二維K-L變換將視網膜圖像三頻帶信息轉換到一個主分量空間，進而降低算法復雜度，并有效地解決了人工通過加權的方式得到灰度圖像的問題，其數據變換矩陣定義如下：

(1)

式中：αk=(αkR,αkG,αkB)T是圖像p的協方差矩陣的特征向量矩陣；k=(1,2,3)；p=(pR,pG,pB)T表示RGB三通道圖像。其主分量信息需要將數據進行協方差對角化，故定義協方差矩陣C(i,j)為：

(2)

(a) 綠色通道 (b) 第一主分量 (c) CLAHE圖2 K-L變換圖像與對比度增強圖像

1.3 偽影矯正

考慮到視網膜圖像包含光照不均勻與偽影等現象，因此利用改進的雙邊濾波Retinex[13]模型。經CLAHE增強的視網膜圖像Ic，可視為反射率圖像Rz與照明圖像L的兩個分量分塊乘法：

Ic=Rz×L

(3)

本文定義x作為輸入圖像Ic的像素值，反射率圖像Rz(x)的像素值x可通過計算輸入圖像Ic(x)和照明圖像L(x)的對數之差來獲得，Rz(x)定義如下：

Rz(x)=log(Ic(x)+1)-log(L(x)+1)

(4)

L(x)定義為：

(5)

歸一化因子M(x)定義為：

(6)

式中:參數φ表示用于測量輸入圖像像素Ic像素與相鄰像素之間空間關系的窗口大小，設置窗口大小為3×3像素[13]。

空間關系g(,x)定義為：

(7)

相似性關系s(,x)定義為：

(8)

式中：d(·)為歐氏距離；參數σd與σr分別表示視網膜圖像空間擴展與強度擴展。函數g(,x)通過歐式距離測量Ic與空間像素關系；函數s(,x)代表Ic與強度相似性關系。

該模型的雙邊濾波法平滑視網膜圖像，保持小鄰域內像素x的強度值彼此近似，并且濾波強度基本一致，同時歸一化項M確保所有臨近強度的權重之和為1，以此充分保留血管的紋理特征，如圖3所示。

(a) 傳統Retinex (b) 本文Retinex圖3 視網膜圖像Retinex矯正

由圖3可看出：傳統Retinex矯正雖然整體亮度得到提高，并且抹掉大部分視網膜圖像的偽影，但仍然存在大面積灰度信息不均勻的光照影響。而本文Retinex算法較好地矯正了視網膜圖像的偽影現象，同時保留血管紋理特征等結構信息，較好地解決后期血管分割將黃斑和視盤誤分割為血管的現象。

2 特征提取

2.1 形態學濾波

利用多尺度形態學Top-Hot變換[14]，在圖像整體增強的同時，提取視網膜血管的微血管特征信息。通過控制圖像邊緣梯度信息控制因子wi，調整相鄰血管像素尺度的差值，降低視盤與黃斑等特征信息的干擾，提高視網膜圖像微小血管的多尺度亮、暗細節特征。多尺度Top-Hot模型定義如下：

(9)

式中：k是視網膜圖像細節增強因子；Ir是輸入經雙邊濾波的Retinex增強的視網膜圖像；fT是多尺度Top-Hot增強后圖像;Dopi與Dcli分別為視網膜圖像亮細節與暗細節特征?？刂埔蜃觲i的值為：

(10)

式中:x、y為經預處理后的視網膜圖像鄰域像素值；gi max與gi min分別為gi的最大值與最小值；gi是視網膜圖像膨脹與腐蝕之差。wi的變化主要由視網膜圖像的梯度信息決定，經多尺度Top-Hot增強的圖像，再利用不同尺度下形態學低帽[15]變換提取血管特征，共計6維特征向量。

2.2 匹配濾波

多尺度高斯匹配濾波器[16]可以用來模擬眼底視網膜血管的整體走向與灰度信息，解決單一尺度下血管寬度變化較大時，血管匹配不全的弊端，其定義如下：

(11)

尺度函數Pi,j(x,y)由兩個不同的響應函數產生，定義如下：

Pi,j(x,y)=Ri(x,y)Rj(x,y)

(12)

多尺度高斯濾波器的響應函數Ri(x,y)定義如下：

Ri(x,y)=fi(x,y)×im(x,y)

(13)

式中：im(x,y)是歸一化的Retinex圖像，其定義如下:

(14)

在本文的多尺度高斯濾波器中，ts為常數項t和濾波器尺度s的乘積；P1,2(x,y)與P3,4(x,y)別用來提取微小血管和寬的血管，歸一化時設定0.75

2.3 Frangi濾波

基于Hessian矩陣的Frangi濾波器[18]，具有可將視網膜寬度不同的血管突出顯示，并且由二階導數在兩個不同尺度上加以分析。Hessian矩陣特征值之間的差異可用于進一步的血管對比度增強和抑制非血管結構?；贔rangi濾波器的Hessian矩陣F(x)定義如下：

F(x)=maxδf(x,δ)

(15)

式中：感興趣像素區域由x定義；δ是用于計算圖像的高斯導數標準差，即血管的尺度方向；f(·)表示濾波器。原函數與高斯二階偏導的卷積可降低二階導數對圖像噪音的敏感度。Hessian矩陣定義如下:

(16)

式中：Ir為經Retinex后的圖像；Gδ(x,y)為高斯函數；Hessian矩陣中的Hxx、Hxy、Hyx和Hyy表示帶有高斯核函數的二階方向導數。定義λ1與λ2為H(x,y)特征值，由0<|λ1|≤|λ2|確定是否為血管像素和明顯的灰度值變化，同時增強血管內部亮度，矯正血管中心亮線。為此構建血管響應函數Vδ(x,y)以區分背景與血管：

(17)

2.4 2D-Gabor濾波

2D-Gabor[19]濾波使用由正弦平面波調制的高斯核函數，其Gabor核具有不同權重的平行帶，其核參數控制條帶的大小、方向和位置，能夠捕捉到血管空間尺度、位置和相位等局部結構信息，增強血管與背景對比度，且對圖像偽影信息具有良好的適應性。本文給出一種基于小波變換的2D-Gabor濾波，設計如下：

首先，利用帶有變換小波ψb,θ,a的標量積來定義連續小波變換Tψ(b,θ,a)：

(18)

式中：f為視網膜圖像平方可積函數；x為圖像像素值；參數b、θ、a分別表示位移矢量、旋轉角、尺度；r為相位偏移量；ψ*為ψ的復共軛；Cψ與ψ分別為歸一化參數和小波分析。

然后，給出2D-Gabor小波的定義如下：

(19)

(20)

其中:θ在Gabor小波變換的跨度以10°為步長[19]，從0°變化到170°。取每個尺度濾波結果的最大值作為該點的血管特征，該特征能有效區分背景與血管之間的灰度信息，共計12維特征向量。

以上4個濾波器特征圖如圖4所示。由于各個濾波器提取的維數特征不一，故將各個特征向量進行歸一化到[-1,1]，并對每幅DRIVE和STARE的訓練集選取7 500個血管與背景像素點，作為正反樣本點，并將血管像素記為1，背景像素記為0。

(a) 形態學濾波 (b) 匹配濾波 (c) Frangi濾波 (d) Gabor濾波圖4 多尺度濾波特征圖

3 AdaBoost分類器及后處理

3.1 AdaBoost算法步驟

AdaBoost的學習過程是基于定義一個分類器，該分類器最小化訓練集X的預測誤差，這個最小化可以視為在訓練集上沿著一組權重Di=(h1,h2,…,hT)(t=1,2,…,T)，最小化分類錯誤，具體步驟如下：

Stept1將經多尺度濾波器提取的視網膜血管特征向量中血管像素與非血管像素以1∶4的比例因子，作為訓練數據集X=(x1,y1),(x2,y2),…,(xm,ym),yi∈{-1,1}，i=1,2,…,m,m為訓練集數目，設定迭代次數為Tmax。

Stept4更新訓練數據集的權值分布，提升樣本在弱分類器中訓練集錯誤分類的權重。

Stept5經過多次迭代，當所有樣本錯誤率小于50%或者幾乎為0，分類器的精度將不再提升，弱分類器通過檢查所有特征的閾值來搜索強分類器H(x)，最終的強分類器由多個弱分類器加權組合而成，并將閾值應用于分類器的輸出以決定分割的血管與背景像素。H(x)定義如下：

(21)

權重迭代規則定義如下：

(22)

3.2 后處理

AdaBoost分類器在根據數據集特征分割血管時，存在將塊狀小斑點區域誤識別為血管。因此，本文構建一個8連通域的幾何形態算子。該算子的閾值由圖像的灰度統計直方圖確定，通過統計去除小于30個像素組成的非細長和非連通區域，得到最終分割結果。由連通域面積(Area)和連通域寬(W)、高(H)信息構建如下幾何算子[7]：

(2)W×H<3.5Area；

(3)Area<30。

4 AdaBoost分類器及后處理

本實驗仿真平臺為MATLAB R2018a，計算機配置為Intel(R)CoreTMi5- 7300 CPU@2.5 GHz,16.00 GB內存,64-bit Windows 10。

4.1 本文方法血管分割圖像比較

圖5展示了本文算法在健康視網膜圖像與病變視網膜圖像的分割結果。其中圖5(a)與(e)為病變視網膜圖像與健康視網膜圖像灰度圖，(b)與(f)為第二專家手工分割結果。觀察(c)與(g)可知本算法在病變視網膜上會分割出多余的病灶與偽影區域，而在健康視網膜的分割結果比較理想，血管結構清晰，基本沒有多余的類病灶區域的出現。(d)與(h)為由連通域除去小于30像素之后圖像。通過觀察(d)利用后處理的方法，分割結果得到了顯著的改善，從而證明了本文后處理的有效性與必要性。

(a) Image0005 (b) 金標準 (c) 初分割圖 (d) 后處理圖

(e) image0163 (f) 金標準 (g) 初分割圖(h) 后處理圖圖5 本文算法分割結果

4.2 不同文獻血管分割圖像比較

圖6展示的是本文算法、文獻[8]和文獻[20]在DRIVE與STARE數據庫的部分圖像的分割結果。其中：(a)為第二專家金標準圖像；(b)為本文算法分割結果；(c)為文獻[20]利用全連接神經網絡與條件隨機場的血管分割圖像；(d)為文獻[7]融合多特征的隨機森林血管分割算法。

(a) 金標準 (b) 本文算法 (c) 文獻[20] (d) 文獻[8]圖6 不同算法分割結果

由圖6前兩行DRIVE數據集和第三行STARE數據集分割結果可知，文獻[8]和文獻[20]均出現視盤誤分割為血管的現象，同時文獻[20]丟失了大部份微血管信息，文獻[8]血管出現過粗的現象。然而本文算法較好地解決了上述兩種算法存在的缺陷，在視盤處基本不存在誤分割的可能，假陽性較低，從而說明引入偽影矯正和二維K-L變化很好抑制了部分噪音信息，提高算法總體魯棒性。并且本文相比上述算法能夠分割處更多的微血管信息，泛化性較強。觀察第二行分割結果可以看出，本文分割結果較文獻[8]和文獻[20]噪音含量極少，血管紋理結構清晰。

綜上所述，本文算法相比以上算法具有較小的誤分割率，能很好地減少視盤被誤分割可能，在血管交叉處不易斷裂且相鄰血管不易相連，對偽影等病灶處有較好的魯棒性。對于主血管末端的微小血管，本文算法能得到較好的保留，且不易存在微小血管斷裂與分割過粗的現象，有較好的連通性與完整性，在保持精度與特異性高的情況下，盡可能分割出較多的血管。

4.3 視網膜分割評價指標

為了進一步定量分析本文算法對視網膜血管提取的算法性能，定義如下幾個評價指標：

(23)

(24)

(25)

式中：TP、TN、FP、FN分別表示真陽性、真陰性、假陽性、假陰性。敏感度(Sensitivity)又稱真陽性率，表示正確分類血管像素占真實血管像素的百分比；特異性(Specificity)表示正確分類的非血管像素占真實非血管像素的百分比; 準確率(Accuracy)表示正確分類血管和非血管像素占整個圖像總像素的百分率。

4.3.1融合多種濾波器與單類型濾波器的比較

為了驗證本文將AdaBoost分類器與多尺度多類型濾波器結合的有效性，給出在DRIVE數據庫下單一濾波器取得的精度與靈敏度指標，并與本文算法的最終結果進行比較，如表1所示。

表1 各個濾波器與多類型濾波器性能比較

由表1可知單濾波器在AdaBoost下也可以得到較好的性能，但單一濾波器存在弊端，將其融合可以減弱單一濾波器相關缺陷。首先，通過基Hessian矩陣的Frangi可增強血管輪廓，多尺度形態學濾波增強微血管信息，降低粗血管與微小血管的對比度信息；2D-Gabor濾波能捕捉血管空間尺度，獲得其相位信息；多尺度高斯濾波可以整合血管整體走向和灰度信息。然后將初步提取的血管特征信息通過AdaBoost多個弱分類器進一步與金標準比對提取血管特征，并加權組合得到強分類器完成最后血管的分類。該過程將無監督學習與有監督學習進行有機融合，使得總體算法精度從平均94%左右達到96%，靈敏度從平均76%提升到83%，更進一步說明本文算法融合4種多尺度濾波器比單一濾波器在AdaBoost分類器分割結果的有效性。

4.3.2不同視網膜分割算法性能比較

表2和表3給出不同算法與本文視網膜分割算法的性能，其中加粗的數值為分割結果最優異的值。選取DRIVE與STARE數據集中準確率、靈敏度和特異性三個指標作為不同算法間的比較標準。

表2 DRIVE數據集視網膜血管分割結果

表3 STARE數據庫視網膜血管分割結果

由表2可知，現有算法在DRIVE數據集血管分割準確率均低于本文算法，說明本文算法在改數據集中具有較強的分割處血管與背景區域的能力。文獻[7]利用PCA的方法分割血管，其特異性高于本文0.003 1，但其靈敏度與準確率分別低于本文算法0.016 8和0.249，說明文獻[7]雖然將陽性較低，但存在血管分割不全和微血管分割嚴重不足的現象。文獻[8]利用隨機森林的血管分割算法，其靈敏度與本文算法相近，但文獻[8]特異性低于本文算法0.018 8，說明本文算法誤分割率遠低于文獻[8]。

從表3可知，本文算法在準確率和特異性方面高于現有算法。文獻[10]靈敏度僅高于本文算法0.010 3，但另外兩個指標低于本文算法，同時文獻[10]采用的SVM分割算法還面臨核函數選擇和易過擬合的問題。文獻[8]融合多種像素特征并利用隨機森林算法分割血管，其靈敏度與準確率雖與本文算法較相近，但均低于本文算法0.004 9和0.002 4，其中特異性低于本文算法0.011 5，從而說明本文算法能在保證較高準確率與靈敏度的情況下，分割的血管圖像假陽性率更低。

綜上所述，本文算法相比以上算法具有較小的誤分割率，能很好地減少視盤被誤分割可能，在血管交叉處不易斷裂且相鄰血管不易相連，對偽影等病灶處有較好的魯棒性。對于主血管末端的微小血管，本文算法能得到較好的保留，且不易存在微小血管斷裂與分割過粗的現象，有較好的連通性與完整性，在保持精度與特異性高的情況下，盡可能分割出較多的血管，提高其應用價值。

5 結 語

本文利用二維K-L變換提取圖像單通道信息，解決選取通道信息經驗化的問題；利用雙邊濾波Retinex算法矯正圖像的偽影等問題，較好地解決了大部分現存算法對存在偽影嚴重的視網膜圖像分割不足的問題，并且進一步抑制視盤區域的對比度，解決后期被誤分割為血管像素的難題。通過巧妙地選取具有不同屬性的多尺度濾波器與AdaBoost相結合，充分發揮二者的優勢，較好地解決了微血管分割不足的難題，并且對于存在病灶眼底視網膜圖像中也具有較高的魯棒性?？紤]到數據集之間的不同圖像具有不同的特征，本文通過組合這些濾波器，有效地克服單一特征或者濾波器存在的微血管分割不足與血管連通性不佳等問題，還可以進一步使分割方法更加穩健。同時，本文算法仍存在部分微血管分割斷裂的現象，故下一步主要工作能有效重構微血管信息增強算法拓撲結構。