基於動態尺度分配的視網膜血管提取方法及系統與流程
2023-06-16 10:13:11 2
本發明涉及基於動態尺度分配的視網膜血管提取方法及系統。
背景技術：
迄今為止，常用的視網膜血管自動提取算法有：
1、基於視網膜血管追蹤方式的算法。這類方法可以較為完整的提取出視網膜血管的網絡，但是算法複雜度較高，運算量較大。此外，對於一些對比度較低的視網膜血管圖像，這類算法的提取準確度不夠。其中比較典型的視網膜血管追蹤算法是由Tolias在1998年提出的基於模糊C均值聚類算法，是在血管的起始處(視盤)挑選出合適的種子點，並由此對整個視網膜血管網絡進行追蹤。建立起視網膜血管橫截面的一維模型，並通過建立起種子點與血管一維模型之間的模糊相似關係，對所有種子點進行分類，判斷其是否屬於血管，從而完成最終的血管網絡分割提取。此類算法存在的最大弊病是種子點的選取會直接影響提取結果，而且在視盤中選取一個最合適的像素點作為種子點是一個比較繁重的工作。此外，此類算法對於血管的分支點的處理效果欠佳，並且會難以避免地丟失細小血管結構，使得分割結果不夠準確。
2、基於分類器的提取算法。這類方法的主要思想是利用視網膜血管所提供的一些先驗信息，構造出一個合適的分類模型，即通常意義上的分類器，並利用構造好的分類器將視網膜圖像中的像素點進行映射分類。Boyce(1999)介紹了一種監督方式的視網膜血管提取算法。這種基於分類器的識別方式對噪聲極為敏感，所以最終的分類效果不是很好。
3、基於匹配濾波器的提取算法，是最為經典也是使用最為廣泛的一種方式。根據血管的灰度分布特徵，大多數匹配濾波算法均選擇高斯濾波器與視網膜圖像進行匹配濾波，生成的響應具有較高的輸出信噪比。其中，最早使用二維高斯匹配濾波器對視網膜血管結構進行提取的方法是Chaudhuri et al.在1989年提出的。此算法基於一個事實：視網膜血管橫截面的灰度分布是服從高斯分布的。因此，若將建立好的高斯濾波器模板與視網膜血管圖像進行不同方向的匹配濾波，當濾波器的尺度與血管的寬度在一定範圍內相吻合時，相應寬度的血管就得以增強，產生很大的卷積響應。但是文獻中僅採用了一種尺度的匹配濾波器，無法使所有寬度的血管均得以增強，並且不論是血管結構還是非血管結構均對高斯匹配濾波器有很強的響應，這樣會導致錯誤分割像素點的出現。針對經典算法的弊端，多種多尺度匹配濾波方案相繼被提出，這一改進旨在利用不同尺度的濾波器對不同寬度的血管進行增強，並在後續步驟中將其準確提取。Bob Zhang在2010年提出的基於高斯一階導數濾波器的視網膜血管提取算法(MF-FDOG)，是在Chaudhuri et al.所提出算法的基礎上，結合高斯一階導數型匹配濾波器與圖像產生卷積響應的特點，使用了多尺度的匹配濾波器對所有視網膜圖像像素點進行遍歷分類。這種方法雖然在一定程度上克服了臺階邊緣對分割結果造成的影響，但是分割結果中仍有噪聲模式存在，一些被噪聲淹沒的細小血管也未被提取出來。此外，此方法對於血管分支點與連通性的處理也有待提高。Qin Li在2012年提出了基於多尺度濾波響應乘積(MPMF)的方法。此方法利用各尺度響應的乘積將各個寬度的血管信息在尺度域進行混合，不僅對血管進行了有效的增強且抑制了噪聲，同時對分割的血管有較好的寬度估計。但本方法仍然存在一定的局限性：其使用了三個不同尺度的濾波器對視網膜血管進行提取，對於變化較大的血管寬度而言，仍存在某些不能被很好增強的血管寬度，因此難以準確地分割出所有血管；另外對於一些存在病變區域視網膜圖像，還會有錯誤分割像素點的出現。
欲以現有方法得到更高的視網膜血管分割精度，目前難以解決的技術問題是：
1、如何利用多尺度方案使各寬度的血管儘量完整地被提取出來；
2、如何減少錯誤分割點以及抑制噪聲；
3、如何準確地估計血管寬度。
技術實現要素：
本發明的目的就是為了解決上述問題，提供基於動態尺度分配的視網膜血管提取方法及系統，通過對視網膜圖像的預處理、圖像分塊、血管分類、動態尺度分配、多尺度匹配濾波、閾值分割和後處理，實現對視網膜圖像的血管提取，並在剔除複雜非血管結構的同時，避免了對血管寬度的過分估計，實現了更簡單、更準確的視網膜血管提取。
為了實現上述目的，本發明採用如下技術方案：
基於動態尺度分配的視網膜血管提取方法，包括如下步驟：
步驟(1)：視網膜圖像預處理：對彩色視網膜圖像的綠色通道分量進行對比度增強；
步驟(2)：圖像分塊：對預處理後的視網膜圖像分割成設定個數的子圖像；
步驟(3)：血管分類：將每個子圖像中血管分為大中小三類；
步驟(4)：動態尺度分配：動態選擇不同尺度的濾波器對不同寬度血管進行增強；
步驟(5)：多尺度匹配濾波：採用多尺度方案並利用高斯匹配濾波模板和高斯一階導數匹配濾波模板對視網膜血管圖像進行濾波處理；
步驟(6)：閾值處理：提取出血管結構並剔除非血管結構，將所有子圖像的提取結果進行重新拼接，得到視網膜血管網絡二值圖像；
步驟(7)：後處理：對閾值處理後的視網膜血管網絡二值圖像進一步消除噪聲，對血管邊緣進行平滑處理，消除圖像中殘留的視網膜邊界，對不連續的細小血管進行斷點連接，從而保留圖像細節，得到分割精度高的視網膜血管網絡圖像。
進一步的，所述步驟(1)的步驟為：
步驟(1-1)：提取彩色視網膜圖像的綠色通道分量：彩色視網膜圖像含有紅色、綠色和藍色三個通道，僅選擇對比度高、噪聲低的綠色通道作為初始處理對象；
步驟(1-2)：多尺度頂帽變換：利用形狀不變，尺寸等差增大的圓形結構元素，對初始處理對象進行頂帽變換處理，增強初始處理對象的對比度；
步驟(1-3)：基於高斯曲線擬合的直方圖線性拉伸：對經多尺度頂帽變換增強所得圖像的灰度直方圖進行基於高斯曲線擬合的直方圖線性拉伸，得到預處理後的視網膜圖像。
進一步的，所述步驟(2)的步驟為：
對預處理後的視網膜圖像按照尺寸進行均等分塊，分割成設定個數的子圖像。
進一步的，所述步驟(3)的步驟為：
根據血管寬度將血管分為大中小三類，依照不同寬度血管的灰度分布特徵，對子圖像包含的血管種類進行判斷。
進一步的，步驟(3)對步驟(2)的子圖像自上而下、自左向右進行編號；先根據子圖像位置的不同將子圖像分為三類：位於原圖四個角的子圖像、位於原圖中央的子圖像和位於原圖四條邊的子圖像；然後針對每類子圖像，根據不同寬度血管的灰度分布特徵，將每個子圖像的血管進行分類。
進一步的，針對每類子圖像，根據不同寬度血管的灰度分布特徵，將每個子圖像的血管進行分類，步驟為：
(3-1)對於位於原圖四個角的子圖像，子圖像中均僅含有中血管；
(3-2)對於位於原圖中央的子圖像，
判斷是否存在大血管，對位於原圖中央的子圖像灰度直方圖的第一設定灰度分布範圍內像素點出現的頻率進行統計，當統計值大於第一設定閾值時，則代表該子圖像中含有大血管，反之，不含大血管；
判斷是否存在中小血管，對位於原圖中央的子圖像灰度直方圖的第二設定灰度分布範圍內的灰度分布進行高斯曲線擬合，得到高斯函數峰值和方差：
若高斯函數峰值大於第二設定閾值且方差小於第三設定閾值，則子圖像第二設定灰度分布範圍內有小血管但無中血管；
若高斯函數峰值小於第二設定閾值且方差小於第三設定閾值，則子圖像第二設定灰度分布範圍內既無小血管又無中血管；
若高斯函數峰值大於第二設定閾值且方差大於第三設定閾值，則子圖像第二設定灰度分布範圍內既有小血管又有中血管；
若高斯函數峰值小於第二設定閾值且方差大於第三設定閾值，則子圖像第二設定灰度分布範圍內無小血管但有中血管。
(3-3)對於位於原圖四條邊的子圖像，
判斷是否存在大血管，對位於原圖邊緣的子圖像灰度直方圖的第一設定灰度分布範圍內像素點出現的頻率進行統計，當統計值大於第四設定閾值時，則代表該子圖像中含有大血管，反之，不含有大血管；
判斷是否存在中小血管的標準，與(3-2)類似，但閾值設定不同，對位於原圖邊緣的子圖像灰度直方圖的第二設定灰度分布範圍內的灰度分布進行高斯曲線擬合，得到高斯函數峰值和方差：
若高斯函數峰值大於第五設定閾值且方差小於第六設定閾值，則子圖像第二設定灰度分布範圍內有小血管但無中血管；
若高斯函數峰值小於第五設定閾值且方差小於第六設定閾值，則子圖像第二設定灰度分布範圍內既無小血管又無中血管；
若高斯函數峰值大於第五設定閾值且方差大於第六設定閾值，則子圖像第二設定灰度分布範圍內既有小血管又有中血管；
若高斯函數峰值小於第五設定閾值且方差大於第六設定閾值，則子圖像第二設定灰度分布範圍內無小血管但有中血管。
所述步驟(4)的步驟為：
利用濾波器尺度與血管寬度之間的關係，動態選擇不同尺度的濾波器對不同寬度血管進行增強。
進一步的，所述步驟(4)的步驟為：
步驟(4-1)：提取大血管所採用的濾波器尺度參數為2；
步驟(4-2)：提取中血管所採用的濾波器尺度參數為1和0.7，提取小血管所採用的濾波器尺度參數為0.5和0.2。
所述步驟(5)的步驟為：
採用高斯匹配濾波模板和高斯一階導數匹配濾波模板同時對增強後的視網膜圖像進行匹配濾波，分別得到高斯匹配濾波器與視網膜圖像的卷積響應及高斯一階導數匹配濾波器與視網膜圖像的卷積響應；並進一步計算得到高斯一階導數匹配濾波器與視網膜圖像卷積響應的局部均值響應。
所述步驟(5)的高斯匹配濾波模板：
高斯型匹配濾波模板上坐標為(x,y)的一點，其權重表示如下：
其中，σ1表示高斯匹配濾波器的尺度，L代表與濾波器模板瞬時方向相平行的血管長度，m代表高斯濾波器模板中係數的均值，
式中Q表示構建的高斯濾波器模板中所包含的點數，N表示濾波器模板所在的鄰域，表示高斯濾波器模板中的一個離散點。
在高斯型匹配濾波模板中減去m項，實現對視網膜圖像背景的平滑處理。並且，由於視網膜血管在平面中的各個方向上均有分布，故而在匹配濾波過程中，濾波器模板需要不斷的旋轉以檢測不同方向上的血管。
所述步驟(5)的高斯一階導數匹配濾波模板：
高斯一階導數匹配濾波模板上坐標為(x,y)的一點，其權重表示如下：
其中，σ2表示高斯一階導數匹配濾波器的尺度，L代表與濾波器模板瞬時方向相平行的血管長度。
所述步驟(6)的步驟為：
利用高斯一階導數匹配濾波器與視網膜圖像卷積響應的局部均值響應對每個像素點的分割閾值進行調整，得到整個圖像的閾值矩陣；將整個圖像的閾值矩陣與高斯匹配濾波器和視網膜圖像的卷積響應進行逐點比較，提取出血管結構並剔除非血管結構，將所有子圖像的提取結果進行重新拼接，得到視網膜血管網絡二值圖像；
進一步的，所述步驟(6)的步驟為：
定義A為高斯匹配濾波器與視網膜圖像的卷積響應，為高斯一階導數匹配濾波器與視網膜圖像卷積響應的局部均值響應；
將所有子圖像利用其灰度分布特徵，按照血管與背景的灰度對比度的高低，以及是否含具有臺階邊緣的非血管結構進行分類，對各子圖像進行動態的閾值處理；
通過將整個圖像的閾值矩陣和高斯匹配濾波器與視網膜圖像的卷積響應逐點比較，每個像素點均與對應點的閾值進行比較，大於閾值的像素點被分類為血管，否則被分類為非血管；
最後，將各子圖像的不同尺度分割結果利用邏輯「或」結合起來，將所有子圖像的分割結果進行重新拼接，得到閾值處理後的視網膜血管網絡二值圖像。
所述將所有子圖像利用其灰度分布特徵，按照血管與背景的灰度對比度的高低，以及是否含具有臺階邊緣的非血管結構進行分類的步驟為：
步驟(6-1)：利用步驟(3)血管分類過程中得到的高斯曲線方差作為判斷圖像對比度高低的依據；
方差c小於設定值cth，表示子圖像血管與背景的灰度對比度高，血管提取容易；
方差c大於設定值cth，表示子圖像血管與背景的灰度對比度低，血管提取困難。
對低對比度子圖像需進行如下處理：對第二設定灰度分布範圍進行分段處理，以像素點出現頻率最高的灰度值為界，分為高低兩個灰度級範圍。
步驟(6-2)：
對圖像灰度直方圖第三設定灰度分布範圍內像素點出現的頻率和s進行統計，並設立第七設定閾值，頻率大於第七設定閾值時，代表子圖像中含具有臺階邊緣的非血管結構。
若判斷出含具有臺階邊緣的非血管結構，那麼接下來要對其進行定位。進一步對局部均值響應的各點進行遍歷判斷，當中坐標為(x,y)的某一任意點響應幅度大於bth時，認為(x,y)位置附近將會有具有臺階邊緣的非血管結構像素點出現；小於bth時則認為(x,y)位置附近將會有血管像素點出現。提取步驟(3)中所判定的三種不同規格血管時，所設置的bth值也是不同的：
bth＝1小血管
bth＝1.1中血管
bth＝1.2大血管
步驟(6-3)：對各子圖像的分割閾值由多閾值公式確定：
其中，Tk,j為子圖像的閾值矩陣，hi,j為常係數，為高斯濾波響應A的局部均值響應，k為圖像血管與背景的灰度對比度對閾值個數的影響，j為有無臺階邊緣對閾值個數的影響，K(c,cth)與J(s,sth)的取值具體如下：
其中，c是高斯擬合曲線的方差，當c小於閾值cth時,表明對比度高，在整個灰度級範圍內選擇閾值，此時K的取值為1；反之，在對比度低的情況下，分為高低兩個灰度級範圍選擇閾值，此時K的取值為2。s是第三設定灰度分布範圍內像素點出現的頻率和，如果s小於sth,則不存在臺階邊緣，不需對閾值針對血管結構和非血管結構進行區分，此時J的取值為1；當s大於等於sth時,存在臺階邊緣，血管結構和非血管結構要分別設置不同的閾值，此時J的取值為2。
綜上所述，針對不同子圖像，有以下四類閾值處理手段：
第一類：對於對比度低且存在臺階邊緣的子圖像，採用h11、h12、h21、h22四個閾值，其中h11針對高灰度級範圍內的血管結構，h12針對高灰度級範圍內的非血管結構，h21針對低灰度級範圍內的血管結構，h22針對低灰度級範圍內的非血管結構；
第二類：對於對比度低且不存在臺階邊緣的子圖像，採用h11、h21兩個閾值，其中h11針對高灰度級範圍內的血管結構，h21針對低灰度級範圍內的血管結構；
第三類：對於對比度高且存在臺階邊緣的子圖像，採用h11、h12兩個閾值，其中h11針對第二設定灰度分布範圍內的血管結構，h12針對第二設定灰度分布範圍內的非血管結構；
第四類：對比度高且不存在臺階邊緣，採用h11一個閾值，h11針對整個灰度級範圍內的血管結構。
整個圖像的閾值矩陣為最終可以通過公式(7)得到分割後的血管網絡：
其中，vess表示各像素點的最終取值。即將整幅圖像的閾值矩陣與圖像的高斯濾波響應A逐點比較，每個像素點均利用對應點的閾值判斷一次，大於對應點閾值的像素點被分類為血管，否則被分類為非血管。
步驟(6-4)：將各子圖像的不同尺度分割結果利用邏輯「或」操作結合起來，再將所有子圖像的分割結果進行重新拼接，得到閾值處理後的視網膜血管網絡二值圖像。
所述步驟(7)的步驟為：
首先，利用目標視網膜血管與噪點所在區域的幾何特徵差別，先刪除視網膜血管網絡二值圖像中面積小於設定閾值P的連通域，再利用連通域外接矩形的長寬比，進一步的消除圖像中殘留的噪聲；
其次，使用多尺度的高斯濾波器模版來對閾值分割後的血管邊緣進行平滑處理；
接下來，利用DRIVE資料庫中提供的掩膜圖像，找到視網膜邊界像素點的位置，進而消除圖像中殘留的視網膜邊界；
最後，採用閉運算對視網膜血管網絡二值圖像中不連續的細小血管進行斷點連接。
基於動態尺度分配的視網膜血管提取系統，包括：
視網膜圖像預處理模塊：對彩色視網膜圖像的綠色通道分量進行對比度增強；
圖像分塊模塊：對預處理後的視網膜圖像按照尺寸進行均等分塊，分割成設定個數的子圖像；
血管分類模塊：根據血管寬度將血管分為大中小三類。並依照不同寬度血管的灰度分布特徵，對子圖像包含的血管種類進行判斷。
動態尺度分配模塊：利用濾波器尺度與血管寬度之間的關係，動態選擇不同尺度的濾波器對不同寬度血管進行增強；
多尺度匹配濾波模塊：採用高斯匹配濾波模板和高斯一階導數匹配濾波模板同時對增強後的視網膜圖像進行匹配濾波，分別得到高斯匹配濾波器與視網膜圖像的卷積響應及高斯一階導數匹配濾波器與視網膜圖像的卷積響應；並進一步計算得到高斯一階導數匹配濾波器與視網膜圖像卷積響應的局部均值響應；
閾值處理模塊：利用高斯一階導數匹配濾波器與視網膜圖像卷積響應的局部均值響應對每個像素點的分割閾值進行調整，得到整個圖像的閾值矩陣。將整個圖像的閾值矩陣與高斯匹配濾波器和視網膜圖像的卷積響應進行逐點比較，提取出血管結構並剔除非血管結構，將所有子圖像的提取結果進行重新拼接，得到視網膜血管網絡二值圖像；
後處理模塊：對閾值處理後的視網膜血管網絡二值圖像進一步消除噪聲，對血管邊緣進行平滑處理，消除圖像中殘留的視網膜邊界，對不連續的細小血管進行斷點連接，從而保留圖像細節，得到分割精度高的視網膜血管網絡圖像。
本發明的有益效果：
1、如果只利用高斯濾波器與視網膜圖像所得到的卷積響應進行閾值分割，某些非血管結構就會不可避免地被錯誤分割出來，致使最終的分割精度受到較大影響。因此，本發明使用兩個匹配濾波模版同時作用，來減小非血管結構被錯誤分割的可能性。
2、匹配濾波部分是在對基於高斯一階導數濾波器的視網膜血管提取算法(MF-FDOG)與基於多尺度濾波響應乘積的視網膜血管提取算法(MPMF)的基礎上，將整個視網膜圖像進行分割處理，並針對各子圖像的灰度分布特徵動態的為其分配不同的濾波器尺度。
3、在多尺度方面，為了避免對中小血管的種類產生錯誤的判斷，本發明分別為提取中小血管分配了兩個不同的尺度參數。本發明一共提供了五個尺度參數，分別用於對大、中、小血管進行增強，是目前多尺度方案中設置尺度數最多的算法。使用這種動態分配尺度的多尺度濾波算法，在抑制具有臺階邊緣的非血管結構以及保證血管寬度的同時，使得更多被噪聲淹沒的細小血管得以增強，進而在後續處理中將其分割提取出來，提高了算法的性能。
4、在閾值處理方面，本發明針對各子圖像對比度的不同以及是否包含臺階邊緣這幾種情況進行分類處理。這種根據子圖像不同的灰度分布情況動態的進行閾值處理的方法，使得對比度較低區域中的一些細小血管也能被分割提取出來。
5、後處理方面，本發明在消除圖像噪聲以及錯誤分割點的基礎上，較好的保留了圖像的細節，使血管分割精度得到了進一步的提高。
6、本發明算法具有較低的複雜度且易於實現，對推進視網膜疾病的診斷具有重大意義。
附圖說明
圖1為本發明的流程圖；
圖2為本發明的系統架構圖；
圖3(a)為DRIVE資料庫測試集中的第19號彩色視網膜血管圖像；
圖3(b)為提取出的綠色通道圖像；
圖3(c)為多尺度頂帽變換後的圖像；
圖3(d)為灰度拉伸後的圖像。
圖4(a)為DRIVE資料庫測試集中的第07號彩色視網膜血管圖像；
圖4(b)和圖4(c)為圖像分割後的兩個子圖像；
圖5為各子圖像的編號情況；
圖6(a)為位於原圖四角的子圖像之一；
圖6(b)為位於原圖四角的子圖像之二；
圖7(a1)為中央子圖像a；
圖7(a2)為中央子圖像a的灰度分布直方圖；
圖7(b1)為中央子圖像b；
圖7(b2)為中央子圖像b的灰度分布直方圖；
圖7(c1)為中央子圖像c；
圖7(c2)為中央子圖像c的灰度分布直方圖；
圖7(d1)為中央子圖像d；
圖7(d2)為中央子圖像d的灰度分布直方圖；
圖8為利用高斯擬合曲線的峰值a和方差c來對中小血管劃分示意圖；
圖9(a)為DRIVE資料庫測試集中的第19號彩色視網膜血管圖像；
圖9(b)為預處理後的圖像；
圖9(c)為本發明算法提取的血管網絡；
圖9(d)為ground truth血管網絡圖。
具體實施方式
下面結合附圖與實施例對本發明作進一步說明。
如圖1所示，基於動態尺度分配的視網膜血管提取方法，包括如下步驟：
步驟(1)：圖像預處理
進一步的，所述步驟(1)的步驟為：
步驟(1-1)：提取彩色圖像的綠色通道分量：原始彩色視網膜圖像含有紅色、綠色和藍色三個通道，僅選擇對比度高、噪聲低的綠色通道作為初始處理對象；
步驟(1-2)：多尺度頂帽變換：利用形狀不變，尺寸等差增大的圓形結構元素，對初始處理對象進行多尺度頂帽變換處理，增強初始處理對象的對比度；
步驟(1-3)：基於高斯曲線擬合的直方圖線性拉伸：對經多尺度頂帽變換增強所得圖像的灰度直方圖進行基於高斯曲線擬合的直方圖線性拉伸，得到預處理後的視網膜圖像。
通過預處理步驟增強後的視網膜圖像如圖3(a)、圖3(b)、圖3(c)和圖3(d)所示。
步驟(2)：圖像分塊：對預處理後的視網膜圖像按照尺寸進行均等分塊，分割成設定個數的子圖像；
由於視網膜圖像中血管和背景在局部範圍內都比較均勻，所以將整個視網膜圖像均等的分割成20個等尺寸的子圖像。圖4(a)、圖4(b)和圖4(c)為增強後的視網膜圖像及其子圖像。
步驟(3)：根據血管寬度將血管分為大中小三類。並依照不同寬度血管的灰度分布特徵，對子圖像包含的血管種類進行判斷。
將圖像等分為20個無重疊的子圖，並對分塊自上而下、自左向右進行編號，如圖5所示。根據位置的不同將20個子圖像分為三類；針對血管寬度的不同，將血管分為大中小三類，並在每一個子圖中分別判斷其中的血管種類。
第一類：位於原圖四角的子圖像，編號1、4、17、20。少部分屬於這一類的子圖像圖6(a)所示，其中含有對比度較高的血管，但屬於血管的像素所佔比重較小；多數子圖像均如圖6(b)所示，其中所包含的血管對比度極低，甚至沒有血管像素，所以很難通過分析其直方圖的灰度分布特徵對所含血管的種類進行判斷。針對資料庫中所有圖像進行分析，且為了權衡上述兩種情況，本算法統一為其分配中等規格的尺度。
第二類：處於原圖像中央的子圖像，編號為6、7、10、11、14、15。
此類子圖像處於原圖像中央，因而其中不包含視網膜之外的黑色背景；且圖中的大血管較暗，所以屬於大血管的像素點應處於直方圖的低灰度級範圍內。
1)判斷是否存在大血管。對直方圖第一設定灰度分布範圍內像素點出現的頻率進行統計，當統計值大於第一設定閾值時，則代表該子圖像中含有大血管，反之，不含大血管；圖7(a1)和圖7(b1)分別代表了兩種情況：圖7(a1)低灰度級範圍內存在大血管，圖7(b1)低灰度級範圍內不存在大血管。通過對其相應的灰度直方圖圖7(a2)、圖7(b2)進行對比可以發現：含有大血管的子圖像，在直方圖第一設定灰度分布範圍內像素點出現的頻率明顯高於不含有大血管的子圖像。
2)判斷是否存在中小血管。圖7(c1)和圖7(d1)分別代表了兩種情況：圖7(c1)高灰度級範圍內同時存在中、小血管，圖7(d1)高灰度級範圍內僅存在小血管。通過對其相應的灰度直方圖圖7(c2)、圖7(d2)進行對比可以發現：灰色背景以及中小血管的灰度分布範圍基本集中在第二設定灰度分布範圍之間，且其中小血管與灰色背景的灰度分布範圍較為接近，甚至部分重疊。此時，對直方圖高灰度級範圍內的灰度分布進行高斯曲線擬合，得到的高斯函數峰值a和方差c可以作為區分中小血管的標準。若峰值大，則證明此時子圖像的高灰度級範圍內不僅存在灰色背景還存在小血管，使得在二者灰度分布重疊的範圍內像素點出現的概率較大；若峰值小，則證明此時子圖像的高灰度級範圍內僅存在灰色背景；方差大時，證明在高灰度級範圍內像素點在各灰度級上分布較為分散。因為中血管與小血管或灰色背景的灰度分布範圍重疊較小，所以此時可以認為該子圖像的高灰度級範圍內存在中血管；方差小時，則認為該子圖像的高灰度級範圍內不存在中血管。圖7(c2)和圖7(d2)使上述區分標準得以驗證。
根據以上方差與峰值大小不同的取值，對於子圖像中是否含有中小血管的判斷可以分為以下四類，如圖8所示：
若高斯函數峰值大於第二設定閾值且方差小於第三設定閾值，則子圖像第二設定灰度分布範圍內有小血管但無中血管；
若高斯函數峰值小於第二設定閾值且方差小於第三設定閾值，則子圖像第二設定灰度分布範圍內既無小血管又無中血管；
若高斯函數峰值大於第二設定閾值且方差大於第三設定閾值，則子圖像第二設定灰度分布範圍內既有小血管又有中血管；
若高斯函數峰值小於第二設定閾值且方差大於第三設定閾值，則子圖像第二設定灰度分布範圍內無小血管但有中血管。
第三類：處於原圖像邊緣的子圖像，編號為2、3、5、8、9、12、13、16、18、19。
1)同樣，判斷是否存在大血管。對位於原圖邊緣的子圖像灰度直方圖的第一設定灰度分布範圍內像素點出現的頻率進行統計，當統計值大於第四設定閾值時，則代表該子圖像中含有大血管，反之，不含有大血管；但是由於這一類子圖像分布在原圖邊緣，故而其中一定包含黑色背景，且這些黑色背景的像素點同樣出現在直方圖的第一設定灰度分布範圍內。因此，這一類子圖像用於判斷大血管的閾值大於第二類子圖像中用於判斷大血管的閾值。
2)判斷該子圖像中是否包含小血管或中血管，與第二類子圖像的判斷標準相同。
若高斯函數峰值大於第五設定閾值且方差小於第六設定閾值，則子圖像第二設定灰度分布範圍內有小血管但無中血管；
若高斯函數峰值小於第五設定閾值且方差小於第六設定閾值，則子圖像第二設定灰度分布範圍內既無小血管又無中血管；
若高斯函數峰值大於第五設定閾值且方差大於第六設定閾值，則子圖像第二設定灰度分布範圍內既有小血管又有中血管；
若高斯函數峰值小於第五設定閾值且方差大於第六設定閾值，則子圖像第二設定灰度分布範圍內無小血管但有中血管。
步驟(4)：動態尺度分配：利用濾波器尺度與血管寬度之間的關係，動態選擇不同尺度的濾波器對不同寬度血管進行增強；
所述步驟(4)的步驟為：
由於血管橫截面的灰度分布呈高斯分布，故而將高斯曲線向橫軸上投影，得到的投影長度就是血管的寬度，而正態曲線下[-3σ,+3σ]範圍內的面積超過了99％，由此可以得出血管寬度與濾波器尺度之間的關係為：
σ＝d/6
其中，d表示血管寬度，σ表示濾波器的尺度。
a、本發明的仿真實驗中，提取大血管所採用的濾波器尺度參數為2。
b、由各子圖像的灰度直方圖可以看出，中等寬度的血管與小血管在灰度分布範圍上的界限並不明顯。所以在二者灰度分布重疊的範圍內，很可能會對中小血管的種類產生錯誤的判斷。如若因此為其分配不恰當的尺度參數，將導致部分寬度的血管不能被檢測出來。為避免這一情況的發生，本算法分別為提取中小血管分配了兩個不同的尺度參數。這樣，在整個算法中一共使用了五種尺度參數。以此在減小上述影響的同時，也使得更多寬度的血管得到增強。本發明仿真實驗中，為提取中血管所採用的濾波器尺度參數為1和0.7，提取小血管所採用的濾波器尺度參數為0.5和0.2。
步驟(5)：多尺度匹配濾波；採用高斯匹配濾波模板和高斯一階導數匹配濾波模板同時對增強後的視網膜圖像進行匹配濾波，分別得到高斯匹配濾波器及高斯一階導數匹配濾波器與視網膜圖像的卷積響應；並進一步計算得到高斯一階導數匹配濾波器與視網膜圖像卷積響應的局部均值響應；
本發明採用兩個匹配濾波模板，高斯型匹配濾波模板上坐標為(x,y)的一點，其權重表示如下：
其中，σ1表示高斯匹配濾波器的尺度，L代表與濾波器模板瞬時方向相平行的血管長度，m代表高斯濾波器模板中係數的均值，
式中Q表示構建的高斯濾波器模板中所包含的點數，N表示濾波器模板所在的鄰域，表示高斯濾波器模板中的一個離散點。
在高斯型匹配濾波模板中減去m項，可以實現對視網膜圖像背景的平滑處理。並且，由於視網膜血管在平面中的各個方向上均有分布，故而在匹配濾波過程中，濾波器模板需要不斷的旋轉以檢測不同方向上的血管。
高斯一階導數匹配濾波模板上坐標為(x,y)的一點，其權重表示如下：
其中，σ2表示高斯一階導數匹配濾波器的尺度，L代表與濾波器模板瞬時方向相平行的血管長度。
步驟(6)：閾值處理：將所有子圖像利用其灰度分布特徵，按照血管與背景的灰度對比度的高低，以及是否含具有臺階邊緣的非血管結構進行分類，對各子圖像進行動態的閾值處理；通過將整個圖像的閾值矩陣與高斯匹配濾波器和視網膜圖像的卷積響應逐點比較，每個像素點均與對應點的閾值進行比較，大於閾值的像素點被分類為血管，否則被分類為非血管；最後，將所有子圖像的分割結果進行重新拼接，得到閾值處理後的視網膜血管網絡二值圖像；
閾值處理部分利用圖像的灰度分布直方圖，對每個子圖像進行動態的閾值處理。將所有子圖像利用其灰度分布特徵按照對比度的高低及是否含有臺階邊緣進行分類。
6-1)區分高對比度圖像與低對比度圖像：利用步驟(3)血管分類過程中得到的高斯曲線方差作為判斷圖像對比度高低的依據；
方差c小於設定值cth，表示子圖像血管與背景的灰度對比度高，血管提取容易；
方差c大於設定值cth，表示子圖像血管與背景的灰度對比度低，血管提取困難。
對低對比度子圖像需進行如下處理：對第二設定灰度分布範圍進行分段處理，以像素點出現頻率最高的灰度值為界，分為高低兩個灰度級範圍。
6-2)區分是否存在臺階邊緣：
對圖像灰度直方圖第三設定灰度分布範圍內像素點出現的頻率和s進行統計，並設立第七設定閾值，頻率大於第七設定閾值時，代表子圖像中含具有臺階邊緣的非血管結構。
若判斷出含具有臺階邊緣的非血管結構，那麼接下來要對其進行定位。進一步對局部均值響應的各點進行遍歷判斷，當中坐標為(x,y)的某一任意點響應幅度大於bth時，認為(x,y)位置附近將會有具有臺階邊緣的非血管結構像素點出現；小於bth時則認為(x,y)位置附近將會有血管像素點出現。提取步驟(3)中所判定的三種不同規格血管時，所設置的bth值也是不同的：
bth＝1小血管
bth＝1.1中血管
bth＝1.2大血管
6-3)閾值選擇
對各子圖像的分割閾值由多閾值公式確定：
其中，Tk,j為子圖像的閾值矩陣，hi,j為常係數，為高斯濾波響應A的局部均值響應，k為圖像血管與背景的灰度對比度對閾值個數的影響，j為有無臺階邊緣對閾值個數的影響，K(c,cth)與J(s,sth)的取值具體如下：
其中，c是高斯擬合曲線的方差，當c小於閾值cth時,表明對比度高，在整個灰度級範圍內選擇閾值，此時K的取值為1；反之，在對比度低的情況下，分為高低兩個灰度級範圍選擇閾值，此時K的取值為2。s是第三設定灰度分布範圍內像素點出現的頻率和，如果s小於sth,則不存在臺階邊緣，不需對閾值針對血管結構和非血管結構進行區分，此時J的取值為1；當s大於等於sth時,存在臺階邊緣，血管結構和非血管結構要分別設置不同的閾值，此時J的取值為2。
針對以上分析，對各子圖像的閾值處理主要分為下列四類情況：
第一類：對於對比度低且存在臺階邊緣的子圖像，採用h11、h12、h21、h22四個閾值，其中h11針對高灰度級範圍內的血管結構，h12針對高灰度級範圍內的非血管結構，h21針對低灰度級範圍內的血管結構，h22針對低灰度級範圍內的非血管結構；
第二類：對於對比度低且不存在臺階邊緣的子圖像，採用h11、h21兩個閾值，其中h11針對高灰度級範圍內的血管結構，h21針對低灰度級範圍內的血管結構；
第三類：對於對比度高且存在臺階邊緣的子圖像，採用h11、h12兩個閾值，其中h11針對第二設定灰度分布範圍內的血管結構，h12針對第二設定灰度分布範圍內的非血管結構；
第四類：對比度高且不存在臺階邊緣，採用h11一個閾值，h11針對整個灰度級範圍內的血管結構。
整個圖像的閾值矩陣為最終可以通過下式得到分割後的血管網絡：
其中，vess表示各像素點的最終取值。即將整幅圖像的閾值矩陣與圖像的高斯濾波響應A逐點比較，每個像素點均利用對應點的閾值判斷一次，大於此點閾值的像素點被分類為血管，否則被分類為非血管。
最後，將各子圖像的不同尺度分割結果利用邏輯「或」操作結合起來，再將20個子圖像的分割結果進行重新拼接，得到閾值處理後的視網膜血管網絡二值圖像。
該動態閾值選擇標準基於圖像的灰度直方圖，利用灰度直方圖中體現的對比度和具有臺階邊緣的非血管結構信息指導閾值的選擇，與單閾值選擇公式相比，既考慮了低對比度區域細小血管的提取，又考慮了對非血管結構的抑制，實驗證明取得了更優的提取效果。
步驟(7)：後處理：對閾值處理後的視網膜血管網絡二值圖像進一步消除噪聲，對閾值分割後的血管邊緣進行平滑處理，消除圖像中殘留的視網膜邊界，對不連續的細小血管進行斷點連接，從而保留圖像細節，得到分割精度高的視網膜血管網絡圖像。
首先，利用目標視網膜血管與噪點所在區域的幾何特徵差別，先刪除視網膜血管網絡二值圖像中面積小於設定閾值P的連通域，再利用連通域外接矩形的長寬比，進一步的消除圖像中殘留的噪聲；
其次，使用多尺度的高斯濾波器模版來對血管邊緣進行平滑處理；
接下來，利用DRIVE資料庫中提供的掩膜圖像，找到視網膜邊界像素點的位置，進而消除圖像中殘留的視網膜邊界；
最後，採用閉運算對視網膜血管網絡二值圖像中不連續的細小血管進行斷點連接。
本發明利用公共資料庫DRIVE對提出的視網膜血管分割算法進行測試。此資料庫中共有40張視網膜圖像，是通過眼底照相機拍攝獲得的，其中每一張圖像的空間解析度均為565×584。該資料庫的建立者將資料庫中的40張視網膜圖像均等的分成了訓練集和測試集，同時在這兩個集合中均提供了一位專家手動分割的結果。但是對於測試集來說，它還提供了第二位專家的手動分割結果，並且，多數的算法均是採用測試集中的圖像來進行仿真實驗的。
將該資料庫提供的第一位專家的手動分割結果，作為衡量本發明算法性能的標準(ground truth)。為了對不同視網膜血管分割算法的性能進行比較，定義了(1)TPR，(2)FPR，(3)ACC三個指標來對算法進行定量描述。定義如下：
其中ACC為視網膜血管分割平均準確率，TN為實驗結果中，被正確分割的血管像素個數；BN為實驗結果中，被正確分割的視場(FOV：圖像中視網膜邊界內的部分)中的背景像素個數。TPR定義為正確分類血管像素點的比例，Nvp為視網膜圖像視場中，目標血管像素個數。FPR定義為視場內屬於非血管但被分類為血管的像素點的比例，Nuvp為視網膜圖像視場中，背景像素的個數。本發明對資料庫DRIVE測試集中的全部20張圖像進行了分割，取各性能指標的平均值作為本算法的最終結果，並與其他算法進行了性能比較，結果如表1所示。
表1是使用本發明算法對DRIVE資料庫測試集全部20張圖像進行處理，所得平均性能指標與其他算法的比較結果。圖9(a)-圖9(d)是本發明算法的20個提取結果之一。圖9(a)為DRIVE資料庫測試集中的第19號彩色視網膜血管圖像，圖9(b)為預處理後的圖像，圖9(c)為本發明算法提取的血管網絡，圖9(d)為ground truth血管網絡圖。對比圖9(c)和圖9(d)可以看出本發明算法能夠將一些細小血管提取出來，並有效地抑制了非血管結構，具有較好的血管分割性能。
表1.DRIVE資料庫測試集血管提取結果(with FOV)
如圖2所示，基於動態尺度分配的視網膜血管提取系統，包括：
視網膜圖像預處理模塊：對彩色視網膜圖像的綠色通道分量進行對比度增強；
圖像分塊模塊：對預處理後的視網膜圖像按照尺寸進行均等分塊，分割成設定個數的子圖像；
血管分類模塊：根據血管寬度將血管分為大中小三類。並依照不同寬度血管的灰度分布特徵，對子圖像包含的血管種類進行判斷；
動態尺度分配模塊：利用濾波器尺度與血管寬度之間的關係，動態選擇不同尺度的濾波器對不同寬度血管進行增強；
多尺度匹配濾波模塊：採用高斯匹配濾波模板和高斯一階導數匹配濾波模板同時對增強後的視網膜圖像進行匹配濾波，分別得到高斯匹配濾波器與視網膜圖像的卷積響應及高斯一階導數匹配濾波器與視網膜圖像的卷積響應；並進一步計算得到高斯一階導數匹配濾波器與視網膜圖像的局部均值響應；
閾值處理模塊：利用高斯一階導數匹配濾波器與視網膜圖像卷積響應的局部均值響應對每個像素點的分割閾值進行調整，得到整個圖像的閾值矩陣。將整個圖像的閾值矩陣與高斯匹配濾波器和視網膜圖像的卷積響應進行逐點比較，提取出血管結構並剔除非血管結構，將所有子圖像的提取結果進行重新拼接，得到視網膜血管網絡二值圖像；
後處理模塊：對閾值處理後的視網膜血管網絡二值圖像進一步消除噪聲，對血管邊緣進行平滑處理，消除圖像中殘留的視網膜邊界，對不連續的細小血管進行斷點連接，從而保留圖像細節，得到分割精度高的視網膜血管網絡圖像。
上述雖然結合附圖對本發明的具體實施方式進行了描述，但並非對本發明保護範圍的限制，所屬領域技術人員應該明白，在本發明的技術方案的基礎上，本領域技術人員不需要付出創造性勞動即可做出的各種修改或變形仍在本發明的保護範圍以內。