低尺度血管檢測在視網膜血管分割中的應用*

吳鑫鑫，肖志勇，劉 辰

江南大學 物聯網工程學院，江蘇 無錫214122

1 引言

報告顯示，截至2017 年，全球約4.25 億成人患糖尿病，占總人口的8.8%，預計到2045 年，糖尿病患者可能達到6.29 億。2017 年全球糖尿病患者（20～79歲）分布中，中國達到了1.144 億，排在第一位。據評估，全球有多達2.124 億人或所有20～79 歲糖尿病患者的一半不知道已患病[1]。并且已知病史的糖尿病患者中，視網膜病變的患病率高達65.2%[2]。視網膜血管的形態、寬度、角度等信息可以直接用于疾病的診斷和篩選[3]，因此早期的檢測很重要。

傳統的檢測需要醫生手工標注血管，任務量大，需要耗費近兩小時[4]來完成一張圖片的血管分割。為了節約人力，提高效率，需要借助計算機輔助，自動化地分割血管。視網膜血管分割方法主要有兩大類：一類是非監督學習的分割方法；另一類是監督學習的分割方法。

非監督學習方法主要是根據血管的某個或者多個特性來提取血管。主要有基于形態學的方法、基于血管跟蹤的方法、基于匹配濾波的方法和基于形變模型的方法。其中基于形態學的方法[4-6]，主要利用了血管比背景亮度低的特點，通過形態學膨脹腐蝕等操作，能很好地檢測出大部分血管。但是此類方法對高亮度區域（視盤、光斑）較敏感，并且對比度較小的血管以及細小的血管檢測效果不太理想。基于匹配濾波的方法，主要利用血管具有較小的曲率，并且血管的橫截面近似于高斯曲線的特性。最常見的有高斯濾波[5]、Gabor 濾波[7]，操作簡單，但是對噪音的影響較大。Azzopardi 等人[8]設計出了B-COSFIRE濾波器，該濾波器可以對不同方向上的血管主干和末端進行精確檢測，尤其對細小的血管有很好的檢測效果。基于血管跟蹤的方法[9-10]，主要利用了血管連續的特點，找到起始點，沿著血管網絡不斷延伸，直到滿足終止條件停止。自動跟蹤算法的難點在于初始種子點的設定，人為的設定需要增加人工負擔，自動設定需要精確定位到血管內部。好的跟蹤算法不僅能檢測連續血管，還需要能夠檢測連續的血管段，這樣才不會由于血管的局部連續而過早結束。基于幾何形變模型的方法[11-12]，主要是用連續的曲線來描述血管的輪廓。邊界曲線通過能量函數的參數定義，為了使能量函數達到最小值，邊界會在力的作用下，向血管的邊界靠近，使得邊界上的點處于穩定狀態。這種方法對邊界敏感，因此背景中的相似部分以及噪音會對結果產生一定的影響。

監督學習的方法主要是事先標記好血管點和背景點，再通過構造好的模型來學習輸入到輸出間的映射關系，不斷調整模型。監督學習的改進方法主要有兩個方面：一個是基于輸入的特征，輸入的特征的個數以及特征在血管方面的特征表達能力，Aslani等人[13]選用了高達17 個特征作為輸入，包括了Gabor濾波響應特征、頂帽變換、B-COSFIRE 濾波響應等；而謝林培[14]僅用了彩色圖像的3 個通道作為輸入特征。另一種是基于模型，對于不同的模型，特征提取的能力也不一樣，蔡軼珩等人[4]采用了支持向量機（support vector machines，SVM）作為分類模型，SVM 是比較典型的二分類模型，主要原則為間隔最大化，對眼底圖像的血管類和背景類像素點有很好的分類作用；Vengalil等人[15]提出了一種基于卷積神經網絡（convolutional neural networks，CNN）的血管分割方法。該模型由1 個輸入層、8 個卷積層、3 個池化層和1 個輸出層組成。模型的輸入直接是未作任何預處理的彩色眼底圖像，輸出層以內核大小為1×1 的卷積層替換全連接層，直接輸出分割結果。其中，輸出結果為灰度圖像，還需進一步后處理，通過閾值分割，得到最后的血管分割圖。這種以卷積層取代全連接層的網絡也被稱為全卷積神經網絡（fully convolutional networks，FCN），也是圖像語義分割的常用模型，可以模擬人腦的視覺系統處理方式，分割效果甚至能夠超過人眼。姜玉靜[16]在FCN 的基礎上，采用了在邊緣檢測方面效果突出的HED（holistically-nested edge detection）網絡模型[17]，并通過條件隨機場（conditional random field，CRF）進行后處理，得到的血管分割精度高達95.15%。監督學習的方法，在充足的訓練數據和優異的模型下，分割的結果往往表現得比非監督學習的方法更好。

Fig.1 Flow chart of vessel segmentation圖1 血管分割流程圖

本文采用監督學習的方法，僅采用兩個特征作為輸入：一個是彩色視網膜圖像對比度較高的綠色通道；另一個是經過B-COSFIRE 濾波得到的響應結果。網絡模型采用的是LVD（low-scales vessel detection），每個尺寸通過子網絡ADS（asymmetric depthfixed sub-network）提取血管特征，然后將兩個低尺寸下的血管特征得到的單一輸出結果融合原尺寸下的血管特征，得到輸出結果，再二值化，得到最后的分割結果。

2 本文方法

分割血管的流程如圖1 所示，在訓練階段，需要先對訓練圖像進行特征提取，將得到的兩個特征歸一化后，送入LVD 模型中訓練，并通過人工標注圖像進行不斷自我學習。在測試階段，需要將測試圖像進行同樣的特征提取和歸一化操作，然后在訓練好的網絡模型中測試，得到測試結果。最后還需要二值化得到最后的分割結果。

2.1 特征提取

為了選取最典型的特征，本文除了將彩色視網膜眼底圖像中對比度較高的綠色通道作為一個特征外，還使用B-COSFIRE 濾波器得到血管響應的特征。

B-COSFIRE 濾波器是受到簡單細胞神經元對特定結構敏感的啟發，開發出來檢測圖像中帶狀區域。B-COSFIRE 濾波器的輸入需要通過一組線性對齊的高斯差分濾波器（difference of Gaussian，DoG）得到，這種計算模型是模擬丘腦內外側膝狀體背側核細胞（lateral geniculate nucleus，LGN）感受視覺信號變化的過程[18]。圖2（a）給出了一個檢測垂直血管的B-COSFIRE 濾波器模型，其中垂直的白條是人為設計的血管，實心圓代表DoG 濾波器的過濾區域。以血管的中心黑點為模型的中點，將附近多個DoG 感應區域考慮進來后，整個B-COSFIRE 模型的感受野近似虛線圍成的橢圓區域。最后，通過乘法操作將來自一組DoG 濾波器的響應結合起來，得到最終的濾波結果。DoG 濾波器表達式如下：

Fig.2 B-COSFIRE model and categories圖2 B-COSFIRE 模型與類別

其中，外部高斯函數偏差為σ，內部高斯函數標準偏差為0.5σ[19]。對于給定的灰度圖I，濾波后的結果為cσ(x,y)=|I×DoGσ|+，其中|·|+cσ符號表示去負為零，即整流線性單元（rectified linear unit，ReLU）。然后，通過響應結果cσ(x,y)自動構造B-COSFIRE 濾波器。圖2 中顯示了B-COSFIRE 兩種濾波器，其中圖2（b）是對稱的B-COSFIRE 濾波器模型，可以檢測血管主干；圖2（c）是非對稱的B-COSFIRE 濾波器模型，對血管末端有很好的檢測效果。以圖2（b）為例，以血管中心點（標注為1 的點，下同）和其周圍的考察點（在同心圓上DoG 響應最大點，如標注2、3、4、5 的點，下同）構造一個三元組集合，其中σi為外部高斯函數的標準偏差，(ρi,φi)是第i個考察點以血管中心點為原點的極坐標，一共m個考察點。為了提高各個點位置的容錯性，需要對DoG 濾波響應進行模糊移位操作。首先，模糊操作是計算DoG 濾波器的加權閾值響應的最大值，通過將DoG 濾波器的響應乘以高斯函數。偏差σ′的計算為σ′=，其中α、為常數，ρi為到血管中心點的距離。移位操作是將每個DoG 模糊響應朝著φi相反的位置移動ρi距離，以到達B-COSFIRE 濾波器的中心。移位量為(?xi,?yi)，其中?xi=-ρicosφi，?yi=-ρisinφi。一個三元組(σi,ρi,φi)模糊移位的響應可以表示為：

其中，-3σ′≤x′,y′≤3σ′。最后，求取集合S中所有三元組的模糊移位的幾何平均，得到最終的濾波結果：

為了檢測不同方向的血管，還需要設計不同方向上的B-COSFIRE 濾波器。而B-COSFIRE 濾波器能很容易調整，直接改變角度φi，就可以構造任意方向上的濾波器：

2.2 歸一化

對數據的歸一化操作，不僅可以降低運算量，而且會更快地收斂。本文將原始數據規范到[-0.5，0.5]的范圍，大致分為兩步。第一步，去掉異常值。本文通過3 倍標準差法去掉異常值，先將特征矩陣減去均值，得到新的特征矩陣，然后對特征矩陣中超過3 倍標準差的部分進行調整，使得新的特征矩陣與均值的差的絕對值不超過標準差的3 倍。第二步，縮放。在第一步操作基礎上，對新得到的特征矩陣進行調整，采用最大最小規范化將數據映射到[0，1]的區間上，再減去0.5 得到。由于輸入數據包含了多個特征，需要對每個特征維度進行歸一化，然后再整合，一起作為LVD 的輸入。

2.3 LVD

特征矩陣進行調整，采用最大最小規范化將數據映射到[0，1]的區間上，再減去0.5 得到。由于輸入數據包含了多個特征，需要對每個特征維度進行歸一化，然后再整合，一起作為LVD 的輸入。FCN 是由Long 等人[20]提出來實現端到端的學習，并且輸出和輸入的數據尺寸相同，適用不同尺寸下的圖像在同一網絡模型中獲取到特定的特征。在基于FCN 的網絡架構中，Xie 等人[17]提出了HED 模型，如圖3 左側所示，該模型是在多尺度下，將每個尺度下的單一輸出融合得到最后的結果，主要用來完成圖像邊緣的提取，具有很好的效果。本文簡化了HED 模型，提出了LVD 模型，模型圖顯示在圖3 中間。相比HED 模型，LVD 模型少了兩個尺度，只保留了3 個尺度，分別代表原尺度和兩個低尺度。并且在兩個低尺度下都有一個單一輸出層（single-output）。為了保留原尺寸的細節特征，需要將兩個低尺寸下的single-output 融合原尺度下提取的feature maps，這樣可以在一定程度上保留所有的細節，避免細節丟失。對于每個尺度下的特征提取，子網絡的選取不僅要考慮到卷積核的尺寸，還需要考慮不同的組合。GoogleLeNet設計了Inception模塊，從最初的Incetption v1到Inception v3，卷積核的尺寸越來越小，從Inception v3 到Inception v4，網絡越來越深[21-24]。本文參考了Inception 模型簡化參數和增加網絡深度的思想，設計了非對稱固定深度子網絡（asymmetric depth-fixed subnet，ADS）。該網絡有10 層卷積層，每一層的卷積核都是非對稱的小卷積核，并且深度固定，如圖4 所示。

Fig.3 HED model and LVD model圖3 HED 模型和LVD 模型

Fig.4 ADS structure圖4 ADS 結構

其中，網絡的第一層采用1×2 的卷積核，最后一層采用3×1 的卷積核。對于子網絡ADS-N，則子網絡ADS 中所有卷積核的深度都為N。將子網絡ADS-N 組合到LVD 模型中，可以得到LVD 完整的網絡結構，如圖5 所示。網絡中有3 個特征提取模塊，主要是以尺度來劃分，分為了原尺寸（scale0）模塊、1/2 尺寸（scale1）模塊和1/4 尺寸（scale2）模塊。在原尺寸下，引入一個1×1 的卷積層，對特征進行降維，得到的feature maps 為M。在兩個低尺寸下，用反卷積層（deconv）上采樣，得到兩個尺寸為scale0 的singleoutput。最后將原尺寸下的M個feature maps 與兩個低尺度下的single-output 通過一個連接（concat）層融合，再用一個1×1 的卷積層卷積，得到一張與原圖尺寸相同的二維血管灰度圖，還需要對灰度圖進一步閾值分割，得到最終的血管分割結果。在LVD 中，每個卷積層后跟一個BN（batch normalization）層來加速收斂和一個ReLU 層進行激活。網絡中每一層的參數顯示在表1 中，其中0 ＜M＜N。

Fig.5 LVD structure圖5 LVD 結構

Table 1 LVD parameters表1 LVD 參數

訓練過程中，整個網絡的損失函數采用均方誤差，對于輸出結果y與人工標注y′的均方誤差計算如下：

其中，n代表像素點的個數，代表平方和。為了防止過擬合，在損失函數中引入了L2 正則化項，最終的損失函數Loss定義如下：

其中，w表示網絡中邊上的權重；k表示正則化項的權重。

3 實驗與分析

3.1 數據集

DRIVE[25]是彩色眼底圖庫，一共40 幅圖像，其中7 幅帶有早期糖尿病視網膜病變，其余33 幅均為正常眼底圖像。每幅圖像的像素565×584，分成訓練集和測試集，每個子集有20 幅圖像，每幅圖像對應兩位專家手動分割的結果。本文選取第一位專家的分割結果作為評判標準，以第二位專家的分割結果作為人眼觀察的結果。

數據集來自DRIVE 數據庫，并通過旋轉、縮放和鏡像操作進行擴充。對于旋轉，每張彩色眼底圖像一次旋轉3°，可以生成120 張不同角度的彩色眼底圖像；對于縮放，每張彩色眼底圖像從原尺寸的0.75 倍到原尺寸的1.25 倍，一次增加0.05 倍，可以產生11 個不同尺寸的彩色眼底圖像；對于鏡像，一張彩色眼底圖像可以生成一個鏡像。一張圖像經過上面3個操作，再通過裁剪和填充使圖像尺寸一致后，可以得到2 640張不同的彩色眼底圖像，對于DRIVE 的20張訓練集，擴充后，整個訓練集擁有高達52 800張彩色眼底圖像。

3.2 評價指標

評價指標主要有敏感性(Se)、特異性(Sp)、準確率(Acc)。計算公式如下：

其中，真陽性（true positive，TP）表示將血管正確分割為正類的像素點數；真陰性（true negative，TN）表示將背景正確分為負類的像素點數；假陽性（false positive，FP）表示將背景錯分為正類的像素點數；假陰性（false negative，FN）表示將血管錯分割為負類的像素點數。另外，本文算法的性能也通過受試者工作特征（receiver operating characteristic，ROC）曲線來評價。ROC 曲線以真陽性率（true positive rate，TPR）為縱坐標，假陽性率（false positive rate，FPR）為橫坐標。AUC（area under ROC curve）面積是ROC 曲線與橫軸之間所夾的面積值，AUC 的值越接近1，表明算法的性能越好。

3.3 參數設置與實驗分析

在特征提取階段，B-COSFIRE的參數來自文獻[8]，對于對稱的B-COSFIRE 濾波器，σ=2.4,σ0=3,α=0.7,ρ={0,2,4,6,8} ；而在非對稱B-COSFIRE 濾波器中，σ=1.9,σ0=2,α=0.1,ρ={0,2,4,6,…,22}；B-COSFIRE濾波器可以檢測12 個方向的血管，即nr=12。

訓練階段，N=50，M=10，整個LVD 模型參數的初始化方式是從截斷的正態分布中輸出隨機值，該正態分布的標準差為0.1；損失函數中的正則化的權值k=0.01，batch_size=5，學習率a=0.1，學習率采用指數衰減，衰減系數為0.99，衰減速度為200，迭代50 000 次，采用隨機梯度下降算法更新參數，并采用滑動平均來控制變量更新的速度，控制速度的衰減系數為0.99。實驗采用Tensorflow框架，在GeForce GTX 1080 Ti的兩個GPU 上并行訓練。

Fig.6 ROC curve圖6 ROC 曲線

Fig.7 Histogram of test result圖7 測試結果的直方圖

Fig.8 Partial test results圖8 部分測試結果

在DRIVE 訓練集上進行測試，得到圖6 的ROC曲線以及AUC 值。從圖中可以看出，ROC 曲線靠近左上角，并且AUC 值高達0.978 2，體現了本文算法的可靠性。網絡模型測試的結果是一張二維的灰度圖，體現的是血管分布的概率圖，還需要進一步二值化得到最終的分割結果。在閾值分割前，先將測試結果采用最大最小歸一化法調整到[0，1]的范圍內。為了分析血管的分布情況，取出一張測試結果的直方圖進行分析，如圖7 所示。可以發現，圖中像素值的分布主要集中在0 和0.9 附近，通過最大類間方差算法可以自動地尋找最佳閾值，并進行閾值分割。最后對得到的網絡模型在DRIVE測試集測試，Se、Sp、Acc平均指標分別達到了0.819 2、0.984 2、0.969 5，部分測試圖像的分割結果顯示在圖8 中。其中，第三列為測試結果，第四列為二值化分割結果，最后一列表示的是二值化結果與專家分割圖的對比，其中黃色區域為相交的區域，而紅色區域是本文算法將背景錯誤劃分為血管的部分，綠色部分為未被檢測出的血管部分。以圖8 中的第一張圖為例，截取了對比結果的部分區域進行分析，為了檢測本文算法對血管末端的檢測效果，放大了局部血管末端，如圖9（a）所示，可以看出檢測的管末端基本跟專家分割的結果保持一致，并且血管的曲率、寬度、長度等屬性幾乎一致，也沒有出現離散的點，具有很好的連續性。本文還放大了密集血管附近的區域，如圖9（b）所示，可以看出，對于空間相鄰的血管，并沒有產生血管粘黏的現象，能很好地將血管間隙中的背景分割開來，達到了很好的檢測效果。本文在實驗中還發現，本文算法檢測出了專家標記的血管部分，如圖10 所示。將原圖的綠色通道局部放大，可以看出原圖中靠近視場邊緣有一個細小的血管，如圖10（a）中的白色箭頭指向的區域，該血管對比度較小，肉眼很難看出，而本文算法卻能夠對對比度極低的細小血管進行檢測，體現了本文算法的健壯性。本文算法也有不足的地方，如圖11 所示。取了兩塊地方，一個是血管末稍處，還有一個是視盤區域，如圖11（a）。其中血管末梢很多未被分割出來，而在測試結果中已經分割出來，如圖11（b）所示。這跟本文的二值化方法有關，最大類間方差考慮的是全局，找到一個具體的閾值整體分割，缺少了血管的局部信息，導致二值化后未被正確分割出來，還需要進一步優化二值化方法，更多地考慮血管的局部信息。在對視盤區域分析時發現，將部分背景區域錯誤檢測為血管。因為視盤亮度較高，邊緣亮度較低，導致邊緣區域對比度較高，被誤檢為血管，模型還需要進一步優化，更多考慮對視盤區域的特征學習。

Fig.9 Local magnification圖9 局部放大圖

Fig.10 Unmarked blood vessels detected圖10 檢測出未被標注的血管

Fig.11 Undetected vessels圖11 未檢測出的血管

本文還與其他算法在DRIVE 數據庫上的表現進行對比，結果展示在表2 中。其中，人眼觀察的結果來自于DRIVE數據庫第2位專家的分割結果。文獻[8]僅采用B-COSFIRE 濾波器對血管進行特征提取，在不添加任何分類器的情況下，只靠閾值分割血管，效果接近人眼觀察的結果；而單獨地使用B-COSFIRE濾波特征，并沒有考慮噪音問題，眼底圖像中存在類似血管的背景噪音，單純的閾值分割很難區別噪音，因此在文獻[26]中，加入了RF（random forest）模型作為分類器，對包含B-COSFIRE 濾波特征的5 個特征對象進行訓練，得到了0.960 6 的分割精確度；在本文算法中，利用全卷積神經網絡的原理設計了新型的LVD 模型作為分類器，同樣采用對血管響應較好的B-COSFIRE 響應特征作為訓練對象，分割精確度達到了0.969 5。而在其他方法[17,27-28]中，訓練數據都未采用眼底圖像的B-COSFIRE 濾波特征，HED-CRF 算法直接采用的是未經處理的彩色眼底圖像，分割精確度僅為0.951 5；FCN 和U-Net算法中都是采用了眼底圖像局部的血管塊，缺乏全局信息，分割精確度不高，分別為0.953 3 和0.953 1。對比發現，B-COSFIRE濾波響應特征更能讓模型學習到血管的特征，并且對整個眼底圖像的特征學習要比基于局部血管塊的特征學習更容易獲取到全局特征，分割效果更好。為了比較分類器的好壞，通過AUC指標進行分析，其中采用全卷積神經網絡模型的AUC指標都超過了0.970 0，本文算法更是達到了0.978 2，而采用RF 模型僅有0.963 2，這是因為神經網絡的學習能力要比傳統的機器學習強，更深層次的網絡更能表達圖像的語義特征。并且本文采用的LVD 模型區別于文獻[26-28]的全卷積網絡模型，沒有采用3×3 的卷積核，全部采用非對稱的小卷積核，模型的表現比U-Net 模型還高。在DRVIE 數據庫測試，Se、Sp、Acc、AUC分別達到了0.819 2、0.984 2、0.969 5、0.978 2，達到了先進水平。

Table 2 Comparison of different algorithms on DRIVE表2 不同算法在DRIVE 數據庫中的對比

4 結束語

本文針對自動分割血管設計了端到端的神經網絡模型LVD，采用了彩色眼底視網膜圖像中對比度較高的綠色通道和對血管有較好響應的B-COSFIRE濾波結果作為訓練的兩個特征，在3 個尺度下，分別用子網絡ADS 提取血管特征，并將融合的結果通過最大類間方差獲取最佳閾值后二值化，得到最終的分割結果。在DRVIE 數據庫中，Se、Sp、Acc、AUC分別為0.819 2、0.984 2、0.969 5、0.978 2，體現了本文算法在血管特征提取上較強的學習能力，充分展現了本文算法在輔助醫療診斷的可行性。