基于融合相位特征的視網膜血管分割算法

李媛媛，蔡軼珩，高旭蓉

(北京工業大學 信息學部，北京 100124)(*通信作者電子郵箱caiyiheng@bjut.com)

0 引言

眼底圖像具有獲取簡便、成像清晰等優勢，可以幫助人們盡早地檢測眼部疾病，及時治療。視網膜血管作為眼底圖像的重要組成部分，是唯一可以通過非入侵的方式觀察到的微小血管，而且一些全身性的疾病，如糖尿病、高血壓等，往往會引起視網膜血管發生病變，所以通過對眼底視網膜血管圖像的分析來進行一些疾病的早期檢測和輔助治療是一種非常好的醫療手段。

鑒于視網膜血管分割的重要性，國內外學者提出了許多有效的視網膜血管分割算法，主要分成非監督(例如，血管追蹤[1]、匹配濾波[2]、基于模型[3]等)和監督[4-9]兩類。其中，監督的方法是將眼科專家手工標記的視網膜血管圖像作為標準樣本，并基于給定的血管特征通過訓練學習優化分類器參數得到最終的分類器，然后將此分類器用于眼底視網膜血管的分割。由于監督的方法依據已有的先驗標記信息來設計，所以基于監督方法的血管分割效果普遍優于基于非監督方法的血管分割效果；但是基于監督的分割方法需要根據提取的圖像特征及標準標記信息構造大量訓練樣本，在訓練學習的過程中逐步優化分類器參數，所以其訓練時間較長，分割效率不高，因此，在基于監督的分割方法中，提取優質的圖像特征用于訓練和分割，不僅可以提高分割的準確率，而且可以顯著減少訓練樣本數和訓練時長，提高算法效率。

Soares等[4]提取眼底視網膜圖像的Gabor變換和像素灰度特征，利用貝葉斯分類器實現血管檢測。其中，不同尺度和不同角度的Gabor變換可以對不同寬度不同走向的血管進行檢測，而且Gabor變換可以濾除噪聲且對光照變化不敏感；但是，該方法只考慮了圖像Gabor特征和灰度特征，對細小血管的檢測不夠準確。吳奎等[5]在Gabor特征的基礎上加入組合線性算子提高了分割的準確率。吳奎[6]提取了眼底視網膜圖像的組合位移濾波響應(Combination Of Shifted FIlter REsponses， COSFIRE)特征，采用貝葉斯模型來進行血管分割。其中，對稱和非對稱COSFIRE濾波器的使用可以檢測出連續血管及血管末端，得到更優質的血管特征。slani等[7]提取眼底視網膜圖像的Hessian矩陣、高帽變換、Gabor變換、條帶選擇組合位移濾波響應(Bar-selective Combination Of Shifted FIlter REsponses, B-COSFIRE)特征以及像素灰度特征用于分類。該方法中作者評估了每個特征在混合特征向量中的重要性，其中B-COSFIRE濾波和Hessian矩陣的預測重要性指標較高。Zhu等[8]提取眼底視網膜圖像的形態學、相位一致性、Hessian矩陣、向量散場以及局部特征用于分類。其中，向量散場特征具有較高的特異性，形態學特征有較高的敏感性。朱承璋等[9]提取眼底視網圖像的不變矩、灰度共生矩陣、高斯拉普拉斯(Laplace Of Gaussian， LOG)算子結合高斯二階導、相位一致性及Hessian矩陣特征用于分類。其中，LOG算子結合高斯二階導特征有較高的敏感性，Hessian特征有較高的特異性，相位一致性特征有較高的準確性；但由于相位一致性特征主要用于血管邊緣檢測，提取出的相位一致性特征中血管中心部分未能檢測出來，所以其敏感性低。

上述基于監督的眼底視網膜血管分割算法大多提取較多的圖像特征，對每個像素構造高維的特征向量用于分類器訓練與分割，耗時長。此外，現有文獻中大多直接使用眼底視網膜圖像的相位一致性特征。相位一致性特征雖對血管邊緣信息有很好的檢測效果，但其對血管中心部位檢測不足，不適用于基于像素分類的視網膜血管分割算法。通過對常用的眼底視網膜圖像特征的分析和研究，本文基于預測重要性指標較高的B-COSFIRE特征、Hessian矩陣特征以及對方向性檢測敏感的Gabor變換特征，針對相位一致性特征對血管中心部位檢測不足的問題，提出了基于融合相位特征的視網膜圖像血管分割算法，該算法克服相位一致性特征對血管中心檢測的不足，且提取較少的特征，實現高效的視網膜血管的分割。

1 基于融合相位特征的視網膜血管分割算法

本文算法是由預處理、特征提取、分類器分類以及后處理4個模塊組成，具體算法流程如圖1所示。預處理后，提取視網膜圖像的Hessian矩陣特征、Gabor特征、B-COSFIRE特征以及相位特征，對圖像中每個像素構造一個4D的特征向量，采用支持向量機(Support Vector Machine, SVM)分類器進行訓練與分類，實現了高效的視網膜血管分割。

圖1 本文算法流程

1.1 預處理

本文通過預處理過程增強眼底視網膜圖像中血管與背景的對比度，便于后續的血管提取。在原始的彩色視網膜圖像中：如圖2(a)，血管與背景部分對比度不高，且存在圖像噪聲，不適合進行后續處理。而各類研究表明，綠色通道分量圖像的對比度高于其他兩個通道分量，且完整地保存了視網膜圖像中的血管結構，故本文選取綠色通道分量圖像進行后續圖像處理；圖2(b)為提取的綠色通道分量的圖像。為進一步增強視網膜圖像的對比度，抑制噪聲，本文采用對圖像增強效果較好的受限制對比度自適應直方圖均衡化(Contrast Limited Adaptive Histogram Equalization, CLAHE)[6]作圖像預處理；圖2(c)為對綠色通道分量進行CLAHE的結果。預處理后，眼底視網膜圖像中血管與背景的對比度明顯增強，更適合后續的處理。

圖2 預處理過程結果

1.2 特征提取

1.2.1 Hessian矩陣

根據眼底視網膜圖像中目標血管呈線性管狀結構的特點，本文采用對線性結構敏感的Hessian矩陣來提取像特征。Hessian矩陣是一個由多元函數二階偏導數構成的方陣，它的特征值和特征向量可以很好地描述視網膜血管類似樹杈形狀的拓撲結構，并且通過濾波函數將眼底視網膜圖像中的噪聲去除。Frangi等[10]最先將Hessian矩陣用于眼底視網膜血管增強，并提出線性濾波函數。

由于視網膜血管的直徑存在變化，單一尺度的Hessian矩陣的血管檢測效果不佳，故將Hessian矩陣的差分運算與高斯核函數結合，通過改變高斯函數的標準差，來獲得不同尺度下的Hessian矩陣。根據高斯函數的卷積性質，空間尺度導數Fab可由輸入圖像I與高斯函數的二階倒數卷積得到：

(1)

通過分析Hessian矩陣的各特征值與其對應的線性目標關系，Frangi等[10]提出在σ尺度下的線性濾波公式：

(2)

1.2.2 Gabor變換

在視網膜血管圖像特征中，Gabor變換具有良好的方向和尺度選擇特性，可以檢測不同方向與不同尺寸的血管，而且Gabor變換對光照變化不敏感，可以克服眼底視網膜圖像光照不均的缺點，因此本文選擇Gabor變換進行圖像特征的提取。Gabor變換即窗函數為高斯函數的加窗傅里葉變換，有時可歸為小波變換的一種。本文中Gabor變換的定義同文獻[4]，即分析小波為Gabor小波的連續小波變換。

2D Gabor變換為：

(3)

(4)

在眼底視網膜圖像血管檢測中，Gabor變換參數設置大多同文獻[1]，ε=4，K0∈[0,3]，旋轉角度θ以10°為步長在10°～170°內旋轉，并取不同方向變換的最大值為尺度a下的輸出響應。膨脹尺度a的選擇應與血管相匹配，包括可能的血管寬度。經過實驗，本文取尺度a∈[2,5]，并將各尺度下的Gabor變換結果進行疊加作為最后輸出。

1.2.3 B-COSFIRE

B-COSFIRE濾波模型由Azzopardi等[11]提出，這種模型模擬復雜感受野細胞的工作機制，對血管類的帶狀結構有很好的檢測效果。在眼底視網膜血管檢測中，可以通過對稱和非對稱的B-COSFIRE濾波器來檢測出連續血管和血管末端，提取更完整的視網膜血管特征。B-COSFIRE濾波響應就是由一組高斯差分(Different of Gaussian, DoG)濾波響應的乘積加權幾何平均得到。

濾波公式為：

(5)

B-COSFIRE濾波器通過自動配置過程來實現其對管狀結構的選擇性，圖3為B-COSFIRE濾波器自動配置模型。在濾波器自動配置過程中，給定一個標準差為σ的DoG濾波器，給定中心點，如圖3(a)、(b)中點1，該濾波器會沿著同心圓進行濾波響應，其中響應最大的位置就是表征強度變化大的點，即感興趣的點，如圖3(a)中的點2、3、4、5，圖3(b)中的點2、3。這些感興趣點的數目取決于所考慮的同心圓的數目和指定的原型模型，并用(σi,ρi,φi)來表示這些感興趣的點，S={σi,ρi,φi|i=1,2,…,n}為這些點的集合，其中：σi為最大的DoG濾波器標準差，(ρi,φi)為這些點相對于中心點的極坐標，n為DoG濾波器個數。

圖3 B-COSFIRE濾波器配置模型

為了提高各個點位置的容錯性，對DoG濾波響應進行模糊和移位操作。模糊操作即計算每一個位置極坐標鄰域的DoG濾波器的最大加權閾值響應，權重為DoG濾波響應與高斯函數Gσ′(x′,y′)的系數乘積，σ′=σ0′+αρi，σ0、α為常數，ρi為與中心濾波器之間的距離。移位操作即將模糊后的DoG濾波響應向φi的相反方向移動距離ρi，模糊移位后的DoG濾波響應為：

yi′)Gσ′(x′,y′)}

(6)

其中：Δxi=-ρicosφi；Δyi=-ρisinφi；-3σ′≤x′，y′≤3σ′。

B-COSFIRE濾波響應公式為：

(7)

1.2.4 相位特征

1)相位一致性特征。

對于一個由方波或三角波組成的傅里葉級數，在方波的階躍或三角波的波峰波谷處，其各次諧波分量的相位相同。故，相位一致性認為圖像信號相位一致性最大處就是圖像特征[12]。在眼底視網圖像血管檢測中，相位一致性特征不受光照強度和圖像對比度的影響，且對于血管的邊緣有很好的檢測效果。

Kovesi[13]提出了二維相位一致性特征的計算公式：

(8)

其中：o表示方向的索引；n為尺度信息；To為噪聲補償；Ano為加權函數；ε為常數，避免分母為零。具體推導過程及參數設置見文獻[13]。

2)融合相位特征。

如圖4(a)所示，相位一致性特征對血管邊緣包括細微的分支血管都有很好的檢測效果，但是無法檢測出血管中心部分。故單一的相位一致性特征不適用于基于像素分類的眼底視網膜血管分割算法。為此，本文提出了一種新的融合相位特征，即將分別提取的相位一致性特征與Hessian矩陣特征進行小波融合，融合后的相位特征不僅保留了良好的邊緣信息，其血管中心部分也被填充為白色，更適用于基于像素分類的視網膜血管的檢測。

在小波融合過程中，源圖像被分解成高頻部分和低頻部分，其中，高頻部分包括源圖像的細節和紋理信息，低頻部分包括源圖像的概貌和平均信息，因此，小波融合更容易提取圖像的結構和細節信息。本文小波融合的具體步驟如下：

步驟1 對兩幅源圖像，如圖4(a)相位一致性特征、圖4(b)Hessian矩陣特征，分別進行二層小波分解，得到各自的高頻子圖像系數與低頻子圖像系數。

步驟2 將兩幅源圖像的高頻系數和低頻系數按照加權平均的融合規則進行融合，得到融合后的高頻系數與低頻系數。

步驟3 對融合后得到的系數進行小波重構，得到融合后圖像，如圖4(c)融合相位特征。融合相位特征即保留了相位一致性特征的良好血管邊緣信息，且血管中間部分被填充為白色。

圖4 相位特征融合過程

本文中提取的4D特征如圖5所示。對視網膜圖像中的每一個像素都構造了4維的特征向量，用于后續的分類。

1.3 SVM分類

本文使用公共的用于血管提取的數字視網膜圖像(Digital Retinal Images for Vessel Extraction, DRIVE)數據庫中的訓練與測試集合，通過LIBSVM工具箱進行SVM訓練和測試。

圖5 本文方法的各特征圖

訓練過程 對DRIVE數據庫的訓練集中的每幅視網膜圖像的每個像素構造4維的特征向量，并分別提取1 500個血管點和非血管點，歸一化后，作為訓練樣本。通過LIBSVM工具箱對SVM進行訓練，得到SVM模型。訓練過程中選擇徑向基函數(Radial Basis Function, RBF)核函數作為核函數。為節省訓練時間，本文中懲罰系數C選擇常用默認值1，gamma選擇0.01，訓練得到SVM模型。

測試過程 從DRIVE數據庫的測試集中任意選擇視網膜圖像并對其像素點提取特征向量，歸一化后，作為測試樣本，輸入訓練后得到的SVM模型，輸出二值化的視網膜血管結構圖。

1.4 后處理

采用連通域度量的方法，將輸出圖像中連通域面積小于25像素的孤立噪聲點去除。

2 實驗結果與分析

2.1 評價指標

為了對算法的分割結果進行評價，將本文的實驗結果與DRIVE數據庫中專家手工分割的結果進行比較，使用準確率(Accuracy, Acc)和受試者特性工作曲線面積(Area Under ROC Curve, AUC)對算法性能進行定量分析。

準確率是評價對視網膜圖像中血管點和非血管點的正確分割的概率，計算公式為：

(9)

其中：TP(True Positive)為真正類，TN(True Negative)為真負類，FP(False Positive)為假正類，FN(False Negative)為假負類。

受試者工作曲線(Receiver Operating Characteristic, ROC)是反映敏感性和特異性連續變量的綜合指標，ROC曲線以真陽性率為縱坐標，假陽性率為橫坐標。受試者特性工作曲線面積(AUC)可由梯形規則計算得到，其值越接近1，表明算法的性能就越好。

2.2 實驗結果

本文以DRIVE數據庫中第1位專家手工分割結果為評價標準，實驗結果如表1所示，本文算法進行血管分割所得平均Acc為0.957 4，平均AUC為0.970 2，結果較好。圖6為本文方法在DRIVE數據庫所有圖像中分割準確率最高的最好分割結果，圖中大量血管均被分割出來，可用于臨床參考。圖7為本文方法在DRIVE數據庫所有圖像中分割準確率最低的最壞分割結果，圖中血管主干部分被分割出來，在發生病變的區域，病變部分被錯誤檢測成血管。雖然圖7中出現成塊的錯誤血管區域，但這依然可以提醒醫生此處異常，為醫生進行疾病檢測節省時間。

表1 本文算法分割結果

圖6 本文方法在DRIVE數據庫中分割準確率最高的最好分割結果

圖7 本文方法在DRIVE數據庫中分割準確率最低的最壞分割結果

為了分析融合相位特征用于像素分類的血管檢測的效果，本文進行了不同特征組合的實驗，結果如表2所示。組合1提取了眼底視網膜圖像的Hessian矩陣、B-COSFIRE濾波和Gabor變換特征，構造3維特征向量用于分類器分割，所得Acc和AUC分別為0.957 3、0.965 8，該特征組合分割效果好。組合2在上述3維特征組合中加入相位一致性特征，構造4維的特征向量用于分類，但血管分割的Acc反而降低。因為相位一致性特征雖對血管邊緣檢測效果好，但其對血管中心部位檢測不足，直接使用相位一致性特征反而降低了分割的準確率。組合3則在上述3維特征組合中加入融合相位特征，即本文方法，實驗得到的血管分割的Acc和AUC均有所提升。結果表明，在組合特征中，融合相位特征比相位一致性特征更適用于基于像素的視網膜血管分割。

表2 不同特征組合分割結果比較

此外，為進一步驗證融合相位特征的有效性，本文將提取的單一的融合相位特征與相位一致性特征，分別通過SVM進行分類，結果如表3所示，單一的相位一致性特征用于像素分類提取出血管的Acc和AUC分別為0.919 1、0.935 9，單一融合相位特征用于像素分類提取出血管的Acc和AUC分別為0.947 8和0.957 8，均高于相位一致性分割結果。實驗結果表明，融合的相位特征比相位一致性特征更適用于基于像素分割的視網膜血管分割算法。

表3 單一特征分割結果比較

為了分析本文算法的有效性，將本文算法與近期基于特征提取且算法準確率較優的文獻算法[6-9,11]，以及近期較熱且具有相似復雜度的基于神經網絡的文獻算法[14-15]進行比較。表4為本文算法與文獻[6-9,11,14-15]算法在DRIVE數據庫上分割準確度的比較，圖8為本文算法與文獻[6,9]算法在DRIVE數據庫上分割結果圖的比較，圖9為本文算法與文獻[7,14]算法在DRIVE數據庫上分割結果圖的比較。如表4所示，與文獻[6-9,11,14-15]算法相比，本文算法分割所得到的平均Acc與AUC總體較高。文獻[9]算法的Acc較高，但其AUC低于本文算法，且如圖8所示，相較于文獻[9]分割得到的血管圖像，本文算法分割得到更多的細小血管且血管連續性較好。文獻[14]算法的AUC較高，但其Acc略低，且如圖9所示，相較于文獻[14]，本文分割雖然部分血管末端連續性不夠，但本文算法分割出更多的細小血管。故從總體指標來看，本文算法結果較優。此外，相較于文獻[8]算法構造的39維特征向量、文獻[9]算法構造的23維特征向量，文獻[14]算法采用的完全卷積神經網絡，本文算法僅提取了4維特征通過SVM分類，降低了算法的復雜度，提高了分割效率。

表4 不同算法分割結果比較

3 結語

針對相位一致性特征對血管中心檢測不足的問題，本文提出基于融合相位特征的視網膜血管分割算法。不同于文獻[8-9]算法，本文算法沒有直接使用相位一致性特征，而是將相位一致性特征同Hessian矩陣特征進行小波融合，得到一種新的融合相位特征。該特征既保留了相位一致性特征中良好的血管邊緣信息，又克服了相位一致性特征中血管中心部位與背景相似的缺點。由于基于像素分類的視網膜血管分割算法需要對每一個像素的特征進行識別分類，因此，本文中的融合相位特征比相位一致性特征更適用于基于像素分類的血管分割算法。此外，本文對每一個像素提取Hessian矩陣、Gabor小波、B-COSFIRE、相位特征，構造4維特征向量用于SVM分類，提取較少的特征與樣本，實現快速、高效的眼底視網膜血管的分割；但是，與現有大多的眼底視網膜血管分割算法一樣，本文對于眼底視網膜中毛細血管的分割不足，損失了一些細小血管的信息。在接下來的工作中，如何對細微小的血管進行更高準確度的分割仍需要繼續研究。

圖8 文獻[6,9]算法與本文算法結果比較

圖9 文獻[7,14]算法與本文算法結果比較