任子暉 邱天舒
(中國礦業大學信息與控制工程學院 江蘇 徐州 221000)
眼底微血管結構復雜,其難點大致分為:① 深度信息的丟失。血管造影時將三維深度信息轉化為二維平面圖像,加劇動靜脈的重疊交叉形成更為復雜的網絡結構。② 采集圖像光線照射不均導致對比度存在差異。③ 病灶形狀多變,影響微血管的精確分割。國內外學者提出眾多的分割方案,大致分為有監督機器學習和無監督算法分割。
有監督法分割精度高且國際研究較廣:文獻[1]采用多尺度2D Gabor函數對輸入圖像進行變換,然后將構造的向量空間輸入貝葉斯分類器識別血管特征;文獻[12]對輸入圖像進行像素分類,提取圖像的灰度特性和不變矩特征并組成7維向量并將征向量作為輸入,通過訓練逐步優化神經網絡參數增強對血管的識別精度。基于機器學習的分割方法需要準備人工分割像素數據作為樣本,實際精度取決于樣本精度且訓練時間較長。
無監督法無需樣本輸入,因此運行效率高,實時性強:文獻[3]提出匹配濾波法,將構造的二維卷積核與圖像進行卷積運算,提取各個方向上的響應最大值得到分割圖像;文獻[4]采用匹配濾波對圖像進行預處理,根據圖像局部屬性運用多個閾值探針檢測血管及背景像素。無監督法能夠提取圖像中大部分血管對于局部低比度區域分割精度較低。
血管跟蹤法是另一種無監督血管分割算法:文獻[5]等人采用了基于粒子濾波法,通過選擇初始種子對毛細血管進行跟蹤,實現血管的細分割;文獻[6]等結合圖像局部灰度特征和連通特性,利用概率追蹤法裁定血管的邊界。邊緣跟蹤法能夠提取相對優質的血管網絡結構,但對于交叉及分叉點的識別準確率欠佳。
針對上述無監督與有監督算法存在的缺陷,本文提出基于離散Curvelet增強和遺傳神經網絡血管分類的新方法。算法流程如圖1所示。

圖1 算法流程圖
傳統空間域法以調整灰度像素值為起點,通過直接修改目標范圍內的灰度值或生成濾波模板與圖像進行卷積運算實圖像增強。由于圖像對比度分布不均容易造成血管信息丟失,因此本文采用離散Curvelet算法,將血管特征進行多尺度、多方向的Curvelet變換,更充分地表現血管的各向異性,更好地描述血管邊緣特征。
Curvelet是一種多分辨率分析工具,本質是將曲線的每一小段近似于直線并進行Curvelet變換,從而更充分地表征曲線特征[7]。
直角坐標系下的二維矩陣:f[t1t2](0≤t1,t (1) 式中:φj,l,k[t1t2]為離散曲波波形。 直角坐標系下的Curvelet徑向窗局部為楔形,整體設置為同心環,如圖2所示。 圖2 離散Curvelet頻域分塊 定義同心環頻域窗函數: (2) (3) (4) 引入相同的斜率: (5) 則有: (6) 式中:Sθl為剪切矩陣,因此離散Curvelet變換定義為: (7) 本文采用基于Wrapping規則的離散Curvelet變換,變換步驟如下: (2) 對每個尺度下的系數矩陣計算:塊最大值系數Mc;噪聲標準差σ。 (3) 不同尺度的系數矩陣乘以相應的增強函數f(x),并通過Curvelet逆變換還原獲得增強圖像。 Curvelet算法能較好地逼近奇異性曲線邊緣。對圖像進行多尺度、多方向的Curvelet分解,展現出不同尺度下的頻率信息。通過改進不同頻率下的非線性增益函數,進一步抑制低頻系數或增強高頻系數,實現血管特征增強[8]。圖像分層尺度如下: S=ceil(log2min(m,n)-3) (8) 式中:m、n為圖像尺度。對于550×550的血管算法設置分解6層:1個Coarse層(內層)方向設置1;4個Detail層(中間層)方向設置16、32、32、64;1個Fine層(外層)方向設置1,如圖3所示。 圖3 Curvelet變換尺度系數重構 Coarse層為17×17低頻系數矩陣,反映圖像背景特征,血管曲線特征較低;Detail層為多矩陣構成的中高頻系數矩陣,包含了血管曲線方向特征,其中第四層的曲線方向特征最明顯;Fine層為550×550高頻系數矩陣,包含了血管輪廓和高頻噪聲信息。 在頻域中對于低頻矩陣系數的增強采用閾值法,即增強大于閾值的低頻系數,抑制低于閾值的頻域系數。低頻增益函數設定為: G(x)=k[g(c(x-b))-g(-c(x+b))] (9) 中高頻及高頻系數中血管輪廓系數和噪聲系數摻雜在一起,線性增強血管系數的同時也增強了噪聲系數。為了克服這一不足,算法將噪聲標準差σ作為分段值,非線性增益函數設定為: (10) 式中:x為中高頻變換系數;k1、k2、k3為自添加的分段函數權重值,通過改變其系數更好地控制分段曲線變換程度,如圖4所示;c′、m為函數分段點調整參數,考慮到噪聲對高頻系數影響,算法設置c′=aσ、m=Kmc。 圖4 中高頻非線性增強曲線 圖4中虛線表示輸入圖像矩陣的頻域信號,實線表示對頻域信號的非線性增強處理,曲線斜率呈下降趨勢,增強高頻系數可減少噪聲干擾。通過閾值法抑制低頻系數實現圖像背景的淡化處理,增強中高頻系數增強血管邊緣信息并抑制噪聲干擾,綜合兩者實現圖像的特征增強。 選擇重建Curvelet中高頻增益函數式(10)的3個分式系數作為變量,通過設定不同的系數來對比不同的增強圖像,以此來獲得最佳參數。參數設置如表1所示。 表1 高頻矩陣變量選擇 圖5給出了對重建的增益函數選擇不同參數產生的增強效果對比圖。 圖5 增強圖像 由圖5可知,不同的參數選擇增強結果也有所不同: (1) (a)、(b)對比:σ主要用來調節圖像清晰度,(a)的噪聲得到抑制但細血管信息也較(b)模糊,(b)的背景噪聲得到放大,所以σ不能太小。 (2) (a)、(c)對比:隨著k1的增加,圖像背景噪聲更加嚴重。因為k1所在的第一段函數斜率較高,k1增加的同時小系數也隨之放大nk1(n>1)倍,k1越大背景噪聲越明顯,所以k1不能過大。 (3) (b)、(d)、(e)、(f)對比:(f)相比于(d)噪聲含量減少,相比于(b)血管更加清晰。原因在于k2所在的第二段函數斜率較k1更為平緩,系數的增加反而能抑制噪聲的影響。綜上,k1選擇1.5,k2選擇2.5,σ選擇0.6。 基于神經網絡的圖像分割需要人工標記的像素數據訓練分類器參數,根據機器學習規則對特征點進行訓練獲得某個特征的像素點屬于血管或噪聲的概率。同時增加多種噪聲樣本輸入提升網絡的容錯能力,最后利用模型對圖像的像素點進行識別輸出血管分類項。 BP神經網絡采用誤差梯度下降學習規則,通過信號的雙向傳輸不斷修正參數進而完成收斂。反復修正參數實現大數據分類是一種耗時的一維空間計算方法,容易造成局部收斂或收斂時間過長。因此需要依靠遺傳算法(genetic algorithm,GA)的全局搜索與并行操作能力,更為精確快速地找到最優解,步驟如下: (1) 人工標記的數據作為樣本輸入,對初始化的權值、閾值進行GA編碼,確定遺傳參數。 (2) 對生成的網絡計算誤差、染色體的適應度值。 (3) 根據改進的選擇算子規則,選擇合適的個體并進行種群遺傳操作,生成新種群。 (4) 若最優適應度值無變化或完成進化上限則停止遺傳操作。 (4) 解碼最優染色體獲得權值、閾值并賦予神經網絡,將樣本數據輸入更新的網絡計算誤差。 (5) 網絡根據誤差梯度下降學習規則進一步更新數據,當誤差達到指定精度或訓練上限則完成訓練。 微血管可視為邊緣平行、有限寬度的線段[9],對血管橫截面進行分析,其灰度像素分布特征函數近似于: K(X,Y)=exp(-x2/2σ2) (11) 根據灰度分布特征函數構造多尺度的2D Gaussian卷積核: f(x,y)=A{1-k*exp(-d2/2σ2)} (12) 由于血管方向呈現360°分布,因此卷積核應進行旋轉以此來匹配多方向上的血管。假設p(x,y)是一個矩陣中的離散點,θi是第i個核的方向,則以(0,0)為中心則旋轉矩陣可以表示為: (13) 基于2D Gaussian卷積核的圖像濾波進一步放大了血管特征,有效地濾除了原RGB圖像中的血管無關項,降低了背景噪聲點在神經網絡中被錯誤識別的概率。將圖像與卷積核進行卷積運算,若卷積核尺度過大則無法區分細血管特征,尺度過小則無法有效濾除噪聲像素易出現血管空洞。因此本文采用多尺度結合[10]方式針對不同管徑進行濾波或特征增強處理,如圖6所示。 圖6 2D Gaussian濾波 2.2.1改進選擇算子 GA中運用較廣的選擇算子大致分為比例選擇法和最優個體保存法[11]。比例選擇法是一種隨機數據采樣法,因此會存在一定的選擇誤差;最優個體法往往只注重少量的最優個體,容易造成種群基因單一,陷入局部收斂。針對以上問題,本文在最優個體保存法的基礎上提出一種新的選擇算子[12],從而在避免陷入局部收斂的情況下找到全局最優解,過程如下: (1) 生成樣本種群,計算個體的適應度值,并從大到小排列。 (2) 對排列的個體按照比例分成優、中、及格、劣四個層次,選擇依據是優者多選增加個體適應度,劣者少選增加種群多樣性。 (3) 針對最優個體法存在的缺陷,對于選擇的個體采用故意制造誤差的方法,優等有90%的概率復制雙份,中等的有10%的概率復制雙份,及格的有0%的概率復制雙份,劣等的直接淘汰。 下面以標有隨機適應度和編號的方框為例: (1) 初始化個體數為12的種群,上層為適應度,下層為編號。 327811961124105123456789101112 (2) 從大到小排列,并按比例分。 121110987654321951164371210128優中及格劣 (3) 優等的有90%概率復制兩份,中等的有10%概率復制兩份,及格的不復制,劣等的淘汰。 121110121110987654或121110987987654 2.2.2個體交叉 種群中的兩個個體根據實數編碼規則,按照一定的概率交叉得到另一個子個體。 (14) 式中:mxy代表第x個染色體的第y位;b為范圍是[0,1]的隨機數。 2.2.3個體變異 變異是隨機從種群中選擇個體,改變多位染色體形成新個體的過程,對第m個個體的n位變異操作如下: (15) 遺傳種群設定為16,進化30代,交叉概率0.25,變異概率0.1,神經網絡最大訓練次數2 000。改進的GA優化BP神經網絡(IGA-BP)對濾波圖像進行分割實驗。同時與未優化GA進行對比,個體適應度曲線如圖7所示,未優化GA終止代數位于30代但未穩定,IGA-BP則在18代趨于穩定,表明改進的GA能夠更為快速地找到權值、閾值。 圖7 進化曲線 訓練誤差如圖8所示,當達到目標精度誤差0.000 01時,未改進算法需要204步達到目標,總耗時65.92 s;改進算法只需要20步就可達到目標,總耗時29.69 s。時間的減少一是由于算子的改進減少了GA全局搜索時間;二是GA為網絡找到初始值較精確的權值、閾值,增加網絡的收斂速度。 圖8 訓練誤差 選取15幅圖片分別使用IGA-BP和未改進遺傳神經網絡(GA-BP)進分割,以準確分割為血管的像素占全局像素的比值為參考對象,分別與人工分割標準數據進行對比,誤差曲線如圖9所示,經改進后的神經網絡與人工分割的誤差相對較小,網絡的誤差預測值與標準數據切合較佳且變化幅度穩定。 圖9 誤差分析 分割算法在DRIVE和STARE數據庫上進行實驗。DRIVE包含40幅眼底圖其中血管占總體圖像的比例分別為12.7%、12.3%;STARE含有20幅眼底圖像其中10幅含病變,血管占總體圖像的比例分別為10.4%、14.0%。 將本文算法和文獻[1]的有監督分類器法、文獻[13]的無監督法進行比較,結果如表2所示,可以看出兩種算法的血管網絡分割結構均不完整。其中:文獻[1]局部細血管區域血管信息損失較嚴重,局部血管連通性較低;文獻[13]較好地保存了大部分血管結構,但是圖像中心視盤未完整分割,保留了非血管信息,局部出現血管粘連現象。 表2 分割對比 續表2 本文算法將頻域方向特征和空間尺度特征相結合,借助神經網絡有效地對血管像素進行分類,在保留血管尾端細節的同時避免了血管的粘連現象,分割對比如圖10所示。 圖10 分割對比 為了進一步比較分割性能,本文從準確率(ACC)、真陽性率(TPR)、假陽性率(FPR)、F-measure 四個方面評價,分別從STARE、DRIVE兩個數據庫隨機選取5幅圖進行質量評估,結果如表3所示。 表3 算法分割效果質量評估 % 本文算法在Chaudhuri的匹配濾波算法作出延伸。由表3數據可以看出,算法的分割準確率較高,性能的提升主要來自離散Curvelet增強與遺傳神經網絡圖像分割結合。其中離散Curvelet主要用來增強血管特征,而遺傳神經網絡則用來提取完整的血管結構。 圖11為本文算法和文獻[1]、文獻[13]的5幅圖像F-measure值比較,對于DRIVE數據庫3種算法的F-measure均值分別為0.764、0.753、0.771。本文算法的F-measure值評價高于文獻13(無監督法),部分值較低的原因在于:① 算法對于病變區域存在誤判;② 部分細血管寬度小于2D Gaussian核函數設定值。 圖11 F-measure值對比 如表4所示,與其他監督類算法進行對比,本文算法分割準確率優于大部分無監督法,低于部分有監督算法。但在TPR方面算法結合了改進的遺傳神經網絡使得在區分噪聲和血管方面有著較優的效果。算法結合了有監督與無監督算法的優點,在保證準確率的前提下有效地降低了對假血管的錯誤分割。 表4 多種算法在DRIVE數據庫的對比 % 本文針對Chaudhuri的匹配濾波法血管結構分割不完整、血管分割模糊的問題提出一種無監督與有監督相結合分割新方法。算法以離散Curvelet和遺傳神經網絡為基礎,首先重建Curvelet增益函數并對系數進行更為靈活的改進,獲得對比度和信噪比更好的增強圖像;然后建立2D Gaussian卷積核對圖像進行濾波增強血管特征;最后運用遺傳神經網絡分析增強圖像的灰度像素輸出血管分類項實現圖像的精確分割。本文算法在提取血管的網絡結構、預防微血管信息丟失方面較佳。血管的精確分割在管徑分析、動靜脈分類、彎曲度分析上具有重要意義,如何進一步抑制邊緣檢測的噪聲并有效地區分假血管提升分割精度是下一步研究方向。


1.2 血管特征Curvelet分解


1.3 系數與增強效果的對比分析



2 GA優化神經網絡的血管分割
2.1 多尺度2D Gaussian圖像濾波

2.2 改進的GA操作



2.3 運行結果分析




3 實驗結果與分析
3.1 效果對比



3.2 分割性能分析



4 結 語