CN109166124B - 一种基于连通区域的视网膜血管形态量化方法 - Google Patents

Abstract

本发明提出一种基于连通区域的视网膜血管形态量化方法，在对眼底图像进行预处理后，获取视网膜血管分割图，然后对血管分割图进行后处理。在此基础上对血管网进行细化和边界化处理，得到血管中心线网络和血管边界图。接着进行角点检测，并将其从血管中心线网络中去除，使血管网中各血管段形成相互分离的连通区域。遍历各血管段，近似血管段中心线，将其变为折线从而计算血管方向。最后计算初始管径值，并在血管段中心线上滑动选取圆心，根据圆心血管方向和前期测得的管径值创建半圆窗口，取该窗口与血管边界的两交点间的距离为新的管径值；由此迭代测得一组管径值，取其中值为该血管段的管径大小。本发明适用于大规模的视网膜血管形态的量化，可靠性高。

Description

一种基于连通区域的视网膜血管形态量化方法

技术领域

本发明属于图像处理技术领域，涉及一种基于连通区域的视网膜血管形态量化方法。

技术背景

随着图像捕获硬件的进步和计算效率的不断发展，加上日益复杂的图像分析和机器学习技术，为获取视网膜等区域生物组织的微小细节提供了基础，使专业人员可通过对图像的检测来发现异常。眼底视网膜血管图像(以下简称眼底图像)是血管系统中唯一可通过无创伤眼底照相直接获取的图像，视网膜中央动脉是全身惟一能够在活体上被直接观察到的小动脉。且相关心血管疾病会影响视网膜血管结构，例如高血压(包括妊娠高血压)、冠心病。因此，观察高血压病患者的眼底情况，常能了解患者机体心、肾、脑等脏器的受害程度，对高血压的诊断及预测有着重要意义。随着近年来高血压、冠心病等心血管疾病的发病率不断的升高，由此引发的视网膜并发症也日益增多。以高血压为例：据美国卫生统计中心报告，高血压影响全球近7000万名美国成年人和大约9亿7000万名成年人。血压升高通常会对身体的血管产生压力，可能导致堵塞或破裂的损害，这些血管的变化也会影响视网膜。高血压视网膜病变是一种常见的心血管系统疾病，临床症状表现为：视网膜动脉管径狭窄变细，管壁反光增强，动静脉管径之比发生变化，从正常的2:3变成1:2～4；动静脉出现交叉压迹现象；血管迂曲；动脉和小动脉分支角度增大；重者动脉呈铜丝状或银丝状。

针对众多的心血管疾病的研究表明，视网膜血管形态变化可能是心血管疾病的重要早期标志。通过观察视网膜血管的各种形态变化，如视网膜动静脉管径、管壁变化、动脉分支角度、交叉压迫现象的有无及程度等，在一定程度上可以了解全身血管的情况，预测冠心病、高血压等心血管疾病的发病率。目前医院里对视网膜血管形态变化的判断主要是根据医学专家的经验，这种方法即耗时耗力又存在主观性的误差，不适合大规模的视网膜血管形态变化的判断。

借助计算机辅助诊断技术，通过开发通用的可重复程序，提供自动和半自动图像分析，对视网膜血管形态进行定量测量(包括血管段识别、血管段的中心线坐标、每个中心线坐标处的局部方向角、血管段管径、血管段的弯曲度以及血管的分支指数等等)，从而对视网膜血管形态变化进行判断，可以协助临床医生提高疾病分析和诊断的可靠性和准确率，且节约医疗资源。另外通过从大量眼底图像中提取的视网膜血管形态定量测量值，可用于发现流行病学间的关联。总而言之，视网膜血管形态的量化对视网膜病变诊断的起到了重要作用，对眼底疾病或者相关心血管疾病的早期发现、辅助诊断有着重要的意义。

发明内容

本发明所解决的技术问题是，针对现有技术的不足，提供一种基于连通区域的视网膜血管形态量化方法，适用于大规模的视网膜血管形态的量化，便于视网膜血管形态的判断，可靠性高。

本发明根据视网膜血管形态的拓扑结构，利用视网膜血管网的交点来遍历确定血管的连通区域，在血管连通域的基础上测量血管段的管径、方向，以解决视网膜血管的量化问题。为了保证血管段连通区域的均衡性，通过Harris角点检测的方式筛选特征点并将其从视网膜血管中心线图中去掉，使得各血管段从特征点处互相分离，形成独立的连通区域，方便下一步对血管的量化。在连通域的基础上，采用道格拉斯-普克算法计算血管段的方向，随后通过创建滑动窗口测量血管段的管径值。使得血管段量化的过程是全自动的，大大提高了量化操作的简便性，且保证了算法在量化视网膜血管时的精准度。

本发明的技术解决方案如下：

一种基于连通区域的视网膜血管形态量化方法，包括以下步骤：

步骤1：输入待量化的眼底图像I_src，对其进行分割，将血管区域像素点的灰度值置为1，非血管区域像素点的灰度值置为0，得到二值化的视网膜血管分割图I_seg，很好的保留了眼底图像中的细小血管；

步骤2：对视网膜血管分割图I_seg进行后处理，去除图中噪声，并填补部分血管中因中心线反射光太强导致分割时产生的空缺，得到血管网状图像I_net；

步骤3：对血管网状图像I_net进行细化处理，得到血管中心线网络；对血管中心线网络在八连通模式下的连通区域进行标记，然后根据各个连通区域内像素的个数来去除较小的连通区域(包含像素点个数小于设定阈值N1的连通区域，得到新的血管中心线网络图I_skl；其中预设阈值N1为经验参数；

步骤4：对血管网状图像I_net进行边界化处理，根据像素的四邻域情况移除血管内部的像素(若一个像素的四邻域像素值都为1，则将该像素置为0)，得到血管边界图I_edge；

步骤5：对I_skl进行角点检测，标记出血管中心线网络中的血管交叉点和分支点，去除此处像素点使血管中心线网络中的血管段相互独立分离，形成不同形状不同长度的连通区域；遍历所有血管段连通区域，确定各连通区域内的像素点个数和坐标信息，进一步去除较小的连通区域(包含像素点个数小于设定阈值N2的连通区域)，得到血管段中心线图I_conn；

步骤6：遍历I_conn中的各血管段，获取血管段的起始点和终止点坐标，采用Douglas-Peucker算法近似该血管段的中心线，获得血管段中心线上各血管像素点的方向；

步骤7：在I_conn中的血管段上，创建迭代半圆窗口测量血管段管径；步骤如下：

7.1)以血管段中心线上的第一个像素点为中心，在垂直于该像素点处的血管方向的方向上向两侧延伸作直线，直到该直线两侧与血管边界图I_edge各存在一个交点；测量两个交点间的距离，作为初始的管径值r；

7.2)以血管段中心线上的第一个像素点为圆心，以初始的管径值r为半径，向该像素点处的血管方向作半圆窗口，测量该半圆窗口与血管边界图I_edge的两个交点之间的距离，作为新的管径值r；

7.3)以血管段中心线上的下一个像素点为圆心(即圆心在该血管段中心线上滑动)，以管径值r为半径，向该像素点处的血管方向作半圆窗口，测量该半圆窗口与血管边界图I_edge的两个交点之间的距离，作为新的管径值r；并返回步骤7.3)进行迭代，直到圆心选至该血管段中心线上最后一个像素点，结束迭代；

步骤8：计算步骤7中得到的各个管径值的中值，作为该血管段的管径值。

进一步地，步骤1中，由于以往方法中将像素点单纯分为血管和非血管两类，对于小血管、血管交叉处、两条血管桥接处、平行血管等特殊位置没进行特殊考虑，导致分割效果差，部分血管信息遗失，而准确提取这些位置血管困难，因此本发明采用多标签分类方法，将眼底图像划分为背景中心区域、背景边缘区域、粗血管中心区域、细血管区域以及粗血管边缘区域这五类区域，将经过专家人工标定的眼底图像(ground truth)作为训练样本，训练基于U-nets卷积神经网络的可信度模型分类器；将待测眼底图像输入预先训练好的可信度模型分类器，得到待测眼底图像中每个像素点位于上述五类区域的预测概率值(可信度)。依据五个预测概率值生成视网膜血管分割图I_seg，生成步骤如下：对于待测眼底图像中每个像素点，首先计算其位于粗血管中心区域、粗血管边缘区域、细血管区域的预测概率值之和，若上述三个预测概率值之和大于预设阈值T，则该像素点位于主血管区域，否则，则不在；然后通过对其位于细血管区域的预测概率值开1/2次方获得其位于细血管区域的预测概率值的增强变换值，计算其位于粗血管中心区域、粗血管边缘区域的预测概率值与其位于细血管区域的预测概率值的增强变换值之和，同样，若三者之和大于预设阈值T，则该像素点位于细小血管的预估区域，否则，则不在；最后从得到的细小血管的预估区域中提取血管骨架获得细小血管精确区域，并将此细小血管精确区域与主血管区域进行合并，作为血管区域，将其中的像素点的灰度值置为1，其它区域作为非血管区域，将其中像素点的灰度值置为0，得到待测眼底图像二值化的视网膜血管分割图I_seg。其中预设阈值T为经验值。

进一步地，步骤2中，通过遍历I_seg中的连通区域，对于较小的连通区域(包含像素点个数小于设定阈值N3的连通区域)，将其中的像素置为0，从而达到去除噪声的目的；部分血管内空缺的填补方法为：先检测去噪后的I_seg中的孔洞(包括血管围成的封闭区域和血管中的空缺)，将孔洞内像素点的灰度值均置为1得到填补孔洞后的图像I_fill。接着将I_fill与去噪后的I_seg做减法运算，检测两图像之间的差异信息从而分离出填补的部分。然后检测填补部分中的连通区域，对于包含像素点个数大于或等于设定阈值N4的连通区域，认为其为血管围成的封闭区域而非血管中的空缺，遍历其中每一个像素点并将其灰度值设为0。此时，填补部分留下来的像素点即为血管中的空缺部分，将其与去噪后的I_seg进行或运算以对血管中的空缺部分进行初步填补，得到初步填补后的图像。最后再对初步填补后的图像进行形态学闭运算先膨胀后腐蚀)，对血管内部的空缺进一步填补，得到最终的血管网状图像I_net。其中设定阈值N3、N4均为经验参数。

进一步地，步骤3中，使用Zhang并行快速细化算法对血管网状图像I_net进行细化处理，得到中心线网络。

进一步地，步骤5中，采用Harris–Stephens算法进行角点检测；算法检测过程如下：首先分别采用x方向梯度算子和y方向梯度算子对血管中心线网络图进行滤波操作得到x方向滤波图像Ix和y方向滤波图像Iy。随后对血管中心线网络图中每一像素点(i,j)计算相关矩阵

其中Ix2(i,j)＝Ix(i,j)²；Iy2(i,j)＝Iy(i,j)²；Ixy(i,j)＝Ix(i,j)·Iy(i,j)，i、j为像素点(i,j)在Ix和Iy上的行列号；为进一步去除噪声影响，对M(i,j)的元素Ix2(i,j)、Iy2(i,j)和Ixy(i,j)进行高斯滤波，得到新的矩阵M(i,j)。接着计算每个像素点的角点响应值R(i,j)＝det(M(i,j))-k·trace(M(i,j))²，其中，det(M(i,j))表示矩阵M(i,j)的行列式，trace(M(i,j))表示矩阵M(i,j)的迹，k为一个经验常数，取值范围为0.04～0.06；遍历血管中心线网络图中所有像素点的角点响应值，得到最大角点响应值，记为R_max；对血管中心线网络图中的任一像素点(i,j)，若其的角点响应值R(i,j)大于其八连通区域内所有像素点的角点响应值(像素点(i,j)的角点响应值R(i,j)在以其为中心大小为3×3的窗口内所有像素点的角点响应值中是最大的，即R(i,j)为局部极大值)且R(i,j)>0.01R_max，则标记该像素点为角点。

进一步地，步骤6中，通过Douglas-Peucker算法近似血管段中心线，将其变为一定程度上保持了原有形状的折线，进而根据折线中的线段斜率确定血管段中心线上各血管像素点的方向。对于任一血管段中心线，具体步骤为：1)找到其上的首尾两点P₁(x₁,y₁)、P₂(x₂,y₂)，并计算连接首尾两点所得的直线

2)根据点到直线的距离公式找到血管段中心线上到直线l距离最大的点P_max(x₃,y₃)，相应的距离为

其中A、B和C是将直线l的方程写为一般式Ax+By+C＝0时对应的系数；3)如果距离D_max＜D_T，则将该直线作为该血管段中心线的近似，该血管段中心线处理完毕，其中D_T是预先给定的阈值，D_T为经验参数；如果D_max≥D_T，则以P_max为分割点，将该血管中心线分为P₁P_max和P₂P_max两段，并分别对两段分别进行步骤1)～3)的处理；当所有分段都处理完毕时，依次连接各个分割点形成的折线，即可以作为该血管段中心线的近似。

进一步地，步骤8中，在血管段中心线上创建半圆窗口测量血管管径。圆心在血管段中心线上滑动，半圆窗口的截断处直径方向与圆心所在位置的血管方向垂直，窗口与血管边界相交于Ai(x₁,y₁)、Bi(x₂,y₂)两点。计算两交点间的距离即此处血管管径

下一个半圆窗口的半径r_i＝d_i。由此迭代创建半圆窗口获得该血管段上的一组管径值，取该组管径值的中值为该血管段的最终管径。

有益效果：

本发明根据视网膜血管形态的拓扑结构，利用视网膜血管网的交点来遍历确定血管的连通区域，在血管连通域的基础上测量血管段的管径、方向，以解决视网膜血管的量化问题。为了保证血管段连通区域的均衡性，通过Harris角点检测的方式筛选特征点并将其从视网膜血管中心线图中去掉，使得各血管段从特征点处互相分离，形成独立的连通区域，方便下一步对血管的量化。在连通域的基础上，采用道格拉斯-普克算法计算血管段的方向，随后通过创建滑动窗口测量血管段的管径值。使得血管段量化的过程是全自动的，大大提高了量化操作的简便性，且保证了算法在量化视网膜血管时的精准度。本发明适用于大规模的视网膜血管形态的量化，便于视网膜血管形态的判断，可靠性高。

附图说明

图1是本发明的流程图。

图2是眼底图像I_src。

图3是视网膜血管分割图I_seg。

图4是对I_seg进行后处理得到的血管网状图像I_net。

图5是根据I_net计算得到的血管中心线网络图I_skl。

图6是根据I_net计算得到的血管边界图I_edge。

图7是角点检测的结果。

图8是根据I_skl计算得到的血管段中心线图I_conn。

图9是Douglas-Peucker算法的示意图。

图10是血管中心线近似的结果。

图11是滑动半圆窗口的示意图。

图12是滑动半圆窗口测量血管段管径的过程。

具体实施方式

如图1所示，本发明实施例提供一种基于连通区域的视网膜血管形态量化方法，包括以下步骤：对眼底图像I_src进行血管分割处理，采用融合多标签分类和深度学习的方法来获取视网膜血管分割图I_seg；接着对血管分割图I_seg进行后处理，去除噪声得到血管网状图像I_net；然后分别对I_net进行细化和边界化操作，对应得到血管中心线网络图I_skl和血管边界图I_edge；接着使用Harris角点检测算法标记I_skl中的血管交叉点和分支点，并将其从I_skl中去除得到血管段中心线图I_conn，进而获取到不同形态且相互分离的血管段连通区域。随后通过Douglas-Peucker算法近似各血管段中心线，从而计算各血管像素点的方向。在采用直线估计的方法确定血管段初始管径后，再结合血管方向以及血管段像素点的坐标信息，迭代作半圆窗口在血管段上滑动测量管径值。具体描述如下：

1.本发明采用多标签分类方法，将眼底图像划分为背景中心区域、背景边缘区域、粗血管中心区域、细血管区域以及粗血管边缘区域这五类区域。将经过专家人工标定的眼底图像(ground truth)作为训练样本，训练基于U-nets卷积神经网络的可信度模型分类器，其输出是输入图像中各像素点位于所述五类区域的预测概率值(可信度)。然后对待测眼底图像I_src，首先对其进行预处理，包括灰度处理、亮度归一化、限制性对比度直方图均衡化、伽马校正，再将其输入预先训练好的可信度模型分类器，得到待测眼底图像中每个像素点位于上述五类区域的预测概率值(可信度)。依据五个预测概率值生成视网膜血管分割图I_seg，生成步骤如下：对于待测眼底图像中每个像素点，首先计算其位于粗血管中心区域、粗血管边缘区域、细血管区域的预测概率值之和，若上述三个预测概率值之和大于预设阈值T(T为经验值，本实施例中取值为0.5。)，则该像素点位于主血管区域，否则，则不在；然后通过对其位于细血管区域的预测概率值开1/2次方获得其位于细血管区域的预测概率值的增强变换值，计算其位于粗血管中心区域、粗血管边缘区域的预测概率值与其位于细血管区域的预测概率值的增强变换值之和，同样，若三者之和大于预设阈值T，则该像素点位于细小血管的预估区域，否则，则不在；最后从得到的细小血管的预估区域中提取血管骨架获得细小血管精确区域，并将此细小血管精确区域与主血管区域进行合并，作为血管区域，将其中的像素点的灰度值置为1，其它区域作为非血管区域，将其中像素点的灰度值置为0，得到待测眼底图像二值化的视网膜血管分割图I_seg。其中预设阈值T为经验值。本发明实施例中通过实验获得其取值为0.5。本发明分割结果既精确分割了粗血管，还能保留细小血管。

2.对视网膜血管分割图I_seg进行去噪处理，通过查找I_seg中的连通区域，在获取每一个连通区域内像素点数量后，对于像素点数量小于设定阈值N3(本实施例中取为100)的连通区域，遍历其中的每个像素点，并将其像素值变为0，去除图中的零散像素点(即像素点数量小于设定阈值N3的连通区域)，避免其对血管段连通区域检测的干扰。另外由于个别眼底图像存在动脉粥样硬化症状，导致分割后得到的I_seg中部分血管内部由于中心线反射光太强而空缺。为对空缺处进行填补，先检测去噪后的I_seg中的孔洞(包括血管围成的封闭区域和血管中的空缺)，将孔洞内像素点的灰度值均置为1得到填补孔洞后的图像I_fill。接着将I_fill与去噪后的I_seg做减法运算，检测两图像之间的差异信息从而分离出填补的部分。然后检测填补部分中的连通区域，对于像素个数大于或等于设定阈值N4(本实施例中取为100)的连通区域，认为其为血管围成的封闭区域而非血管中的空缺，遍历其中每一个像素点并将其灰度值设为0。此时，填补部分留下来的像素点即为血管中的空缺部分，将其与去噪后的I_seg进行或运算以对血管中的空缺部分进行初步填补，得到初步填补后的图像。最后再使用大小为5×5、值全为1的卷积核(矩阵)对初步填补后的图像进行形态学闭运算(先膨胀后腐蚀)，对血管内部的空缺进一步填补，得到最终的血管网状图像I_net。

3.对血管网状图像I_net进行细化处理，使用Zhang并行快速细化算法得到中心线网络：设点P1是待检测点，P2～P9是P1的八邻域像素点，则细化删除条件为：a、2≤N(P1)≤6，N(P1)为P1的八邻域像素点中灰度值非零的像素点个数；b、A(P1)＝1，A(P1)指的是以P2,P3,…,P9,P2为序时这些点的灰度值从0到1变化的次数；c、p2*p4*p6＝0或者p2*p4*p8＝0；d、p4*p6*p8＝0或者p2*p6*p8＝0，其中p2、p4、p6、p8分别为像素点P2、P4、P6、P8的灰度值；如果同时满足以上四个条件则该点可以删除。这样反复迭代，直到所有的点都被遍历为止，结果如图4所示。对血管中心线网络在八连通模式下的连通区域进行标记，然后根据各个连通区域内像素的个数来去除包含像素点个数小于设定阈值N1的连通区域，得到新的血管中心线网络图I_skl；其中预设阈值N1为经验参数；本实施例中取为100；

4.对血管网状图像I_net进行边界化处理。遍历I_net中的像素，若一个像素的四邻域像素值都为1，则将该像素置为0。由此，位于血管边界上的值为1的像素被保留下来，得到血管边界图I_edge。

5.采用Harris–Stephens算法对血管中心线网络I_skl进行角点检测，角点可以描述为以下几种：灰度的梯度局部最大所对应的像素点；两条及两条以上边缘的交点；图像中梯度值和梯度方向的变化速率都很高的点。I_skl中检测的角点是血管中心线的交叉点和分支点。算法检测过程如下：首先分别采用x方向梯度算子[-1,0,1；-1,0,1；-1,0,1]和y方向梯度算子[-1,-1,-1；0,0,0；1,1,1]对I_skl进行滤波操作得到x方向滤波图像Ix和y方向滤波图像Iy。随后对每一像素点计算相关矩阵

其中Ix2(i,j)＝Ix(i,j)²；Iy2(i,j)＝Iy(i,j)²；Ixy(i,j)＝Ix(i,j)·Iy(i,j)，i、j为像素点(i,j)在Ix和Iy上的行列号；为进一步去除噪声影响，对M(i,j)的元素Ix2(i,j)、Iy2(i,j)和Ixy(i,j)进行高斯滤波，得到新的矩阵M(i,j)。接着计算I_skl中每个像素点的角点响应值R(i,j)＝det(M(i,j))-k·trace(M(i,j))²，其中，det(M(i,j))表示矩阵M(i,j)的行列式，trace(M(i,j))表示矩阵M(i,j)的迹，k为一个经验常数，取值范围为0.04～0.06；遍历血管中心线网络图中所有像素点的角点响应值，得到最大角点响应值，记为R_max；对任一像素点(i,j)，若其的角点响应值R(i,j)大于其八连通区域内所有像素点的角点响应值(即R(i,j)为局部极大值)且R(i,j)>0.01R_max，则标记该像素点为角点。最后再根据角点的八邻域情况对初步检测到的角点进行进一步筛查，若某一角点的八邻域中有3或4个像素点的灰度值不为0，则符合血管中心线上交叉点和分叉点的结构，保留该角点，否则去除。具体方法为：设角点P1及其八邻域像素点的灰度值构成的矩阵为：

采用算子

对其进行滤波操作，即将矩阵P与F进行点乘，若滤波结果x_p＝p₉·1+p₂·1+p₃·1+p₄·1+p₅·1+p₆·1+p₆·1+p₈·1+p₁·0为3或4，则代表角点P1的八邻域中有3或4个像素点的灰度值不为0，保留角点P1，否则去除。

在获得最终的角点集合后，将角点的像素值置为1，使血管中心线网络中的血管段相互分离，各血管段形成独立的连通区域。循环遍历所有连通区域，确定各连通域内的像素个数和坐标信息。若某连通区域中的像素个数少于设定阈值N2(本实施例中取为10)个则将其去除，最终得到血管段中心线图I_conn。

6.遍历I_conn中的血管段，获取处于该血管段上像素的坐标信息，采用Douglas-Peucker算法计算各像素点处的血管方向。算法理论是先对曲线进行采样即取点从而近似血管中心线，将其变为一定程度上保持了原有形状的折线，进而根据折线中的线段斜率确定血管方法。算法执行步骤如下：首先，连接该血管中心线段的首尾两点S、E得到直线SE；然后计算中心线段上像素点到直线SE的距离，找到其中的最大距离D_max；若D_max＜D_T，则该直线段作为曲线的近似，该段曲线处理完毕，其中D_T是预先给定的阈值，本实施例中取为1；若D_max≥D_T，则用获得D_max对应的点M将血管中心线分为SM和ME两段，并分别对两段进行前3步的处理；当所有曲线都处理完毕时，依次连接各个分割点形成的折线，即可以作为血管中心线的近似。由此血管中心线段上的每个像素点都在折线上有个近似点，此近似点是过像素点作折线的垂线获得的垂足位置，取近似点所在直线的斜率为该血管像素点处的方向。

7.任意取I_conn中要测量的血管段seg上的一个像素点，以该点为中心，在垂直于该点处血管方向的方向上朝两侧延伸作直线，直到该直线与血管边界图I_edge相交于A(x₁,y₁)、B(x₂,y₂)。计算两交点间的距离

8.取血管段seg上起始点为圆心，设置初始半径r＝d，根据圆心点处的血管方向k_o作半圆窗口，使窗口的半圆截断处直径与血管方向垂直。获取半圆窗口与血管边界图I_edge产生的两交点间的距离d_i，更新半径r_i＝d_i/2+2，圆心在血管中心线上滑动，迭代创建半圆窗口，获得该血管段上的管径值向量。最后取所测得的管径值的中值为该血管段seg的管径值。

Claims (8)

1.一种基于连通区域的视网膜血管形态量化方法，其特征在于，包括以下步骤：

步骤1：输入待量化的眼底图像I_src，对其进行分割，得到视网膜血管分割图I_seg；

步骤2：对视网膜血管分割图I_seg进行后处理，去除图中噪声，并填补部分血管中因中心线反射光太强导致分割时产生的空缺，得到血管网状图像I_net；

步骤3：对血管网状图像I_net进行细化处理，得到血管中心线网络图；对血管中心线网络在八连通模式下的连通区域进行标记，然后去除包含像素点个数小于设定阈值N1的连通区域，得到新的血管中心线网络图I_skl；

步骤4：对血管网状图像I_net进行边界化处理，根据像素点的四邻域情况移除血管内部的像素点，得到血管边界图I_edge；

步骤5：采用Harris-Stephens算法对血管中心线网络图I_skl进行角点检测，标记出其中的血管交叉点和分支点，去除血管交叉点和分支点处像素点，形成多个相互独立的血管段连通区域；遍历所有血管段连通区域，进一步去除包含像素点个数小于设定阈值N2的血管段连通区域，得到血管段中心线图I_conn；

步骤6：遍历I_conn中的各血管段连通区域，通过Douglas-Peucker算法获取血管段中心线和血管段中心线上各像素点处的血管方向；

步骤7：在I_conn中的血管段上，创建迭代半圆窗口测量血管段管径；步骤如下：

7.1)以血管段中心线上的第一个像素点为中心，在垂直于该像素点处的血管方向的方向上向两侧延伸作直线，直到该直线两侧与血管边界图I_edge各存在一个交点；测量两个交点间的距离，作为初始的管径值r；

7.2)以血管段中心线上的第一个像素点为圆心，以初始的管径值r为半径，向该像素点处的血管方向作半圆窗口，测量该半圆窗口与血管边界图I_edge的两个交点之间的距离，作为新的管径值r；

7.3)选取血管段中心线上的下一个像素点为圆心，以管径值r为半径，向该像素点处的血管方向作半圆窗口，测量该半圆窗口与血管边界图I_edge的两个交点之间的距离，作为新的管径值r；并返回步骤7.3)进行迭代，直到圆心选取至该血管段中心线上最后一个像素点，结束迭代；

步骤8：计算步骤7中得到的各个管径值的中值，作为该血管段的管径值。

2.根据权利要求1所述的基于连通区域的视网膜血管形态量化方法，其特征在于，步骤1中，将眼底图像划分为背景中心区域、背景边缘区域、粗血管中心区域、细血管区域以及粗血管边缘区域这五类区域，将经过专家人工标定的眼底图像作为训练样本，训练基于U-nets卷积神经网络的可信度模型分类器；

将待测眼底图像输入训练好的可信度模型分类器，得到待测眼底图像中每个像素点位于上述五类区域的预测概率值；依据五个预测概率值生成视网膜血管分割图I_seg，生成步骤如下：

对于待测眼底图像中每个像素点，首先计算其位于粗血管中心区域、粗血管边缘区域、细血管区域的预测概率值之和，若上述三个预测概率值之和大于预设阈值T，则该像素点位于主血管区域，否则，则不在；然后通过对其位于细血管区域的预测概率值开1/2次方获得其位于细血管区域的预测概率值的增强变换值，计算其位于粗血管中心区域、粗血管边缘区域的预测概率值与其位于细血管区域的预测概率值的增强变换值之和，同样，若三者之和大于预设阈值T，则该像素点位于细小血管的预估区域，否则，则不在；最后从得到的细小血管的预估区域中提取血管骨架获得细小血管精确区域，并将此细小血管精确区域与主血管区域进行合并，作为血管区域，将其中的像素点的灰度值置为1，其它区域作为非血管区域，将其中像素点的灰度值置为0，得到待测眼底图像二值化的视网膜血管分割图I_seg。

3.根据权利要求2所述的基于连通区域的视网膜血管形态量化方法，其特征在于，步骤2中，对视网膜血管分割图I_seg进行后处理，去除图中噪声的方法为：通过遍历I_seg中的连通区域，对于像素点数量小于设定阈值N3的连通区域，将其中的像素置为0。

4.根据权利要求3所述的基于连通区域的视网膜血管形态量化方法，其特征在于，步骤2中，填补部分血管中因中心线反射光太强导致分割时产生的空缺的方法为：先检测去噪后的I_seg中的孔洞，将孔洞内像素点的灰度值均置为1得到填补孔洞后的图像I_fill；接着将I_fill与去噪后的I_seg做减法运算，检测两图像之间的差异信息从而分离出填补的部分；然后检测填补部分中的连通区域，对于像素点个数大于或等于设定阈值N4的连通区域，认为其为血管围成的封闭区域而非血管中的空缺，遍历其中每一个像素点并将其灰度值设为0，此时，填补部分留下来的像素点即为血管中的空缺部分，将其与去噪后的I_seg进行或运算以对血管中的空缺部分进行初步填补，得到初步填补后的图像；最后对初步填补后的图像进行形态学闭运算，对血管内部的空缺进一步填补，得到最终的血管网状图像I_net。

5.根据权利要求1所述的基于连通区域的视网膜血管形态量化方法，其特征在于，步骤3中，使用Zhang并行快速细化算法对血管网状图像I_net进行细化处理，得到中心线网络。

6.根据权利要求1所述的基于连通区域的视网膜血管形态量化方法，其特征在于，步骤5中，采用Harris-Stephens算法进行角点检测；算法检测过程如下：首先分别采用x方向梯度算子和y方向梯度算子对血管中心线网络图进行滤波操作得到x方向滤波图像Ix和y方向滤波图像Iy；随后对血管中心线网络图中每一像素点(i，j)计算相关矩阵

其中Ix2(i，j)＝Ix(i，j)²；Iy2(i，j)＝Iy(i，j)²；Ixy(i，j)＝Ix(i，j)·Iy(i，j)，i、j为像素点(i，j)在Ix和Iy上的行列号；再对M(i，j)的元素Ix2(i，j)、Iy2(i，j)和Ixy(i，j)进行高斯滤波，得到新的矩阵M(i，j)；接着计算每个像素点的角点响应值R(i，j)＝det(M(i，j))-k·trace(M(i，j))²，其中，det(M(i，j))表示矩阵M(i，j)的行列式，trace(M(i，j))表示矩阵M(i，j)的迹，k为一个经验常数，取值范围为0.04～0.06；遍历血管中心线网络图中所有像素点的角点响应值，得到最大角点响应值，记为R_max；对血管中心线网络图中的任一像素点(i，j)，若其的角点响应值R(i，j)大于其八连通区域内所有像素点的角点响应值且R(i，j)＞0.01R_max，则标记该像素点为角点。

7.根据权利要求1所述的基于连通区域的视网膜血管形态量化方法，其特征在于，标记出角点后，根据八邻域情况对角点进行进一步筛查，若某一角点的八邻域中有3或4个像素点的灰度值不为0，则保留该角点，否则去除，获得最终的角点集合。

8.根据权利要求1所述的基于连通区域的视网膜血管形态量化方法，其特征在于，步骤6中，通过Douglas-Peucker算法近似血管段中心线，将其变为折线，进而根据折线中的线段斜率确定血管段中心线上各像素点处的血管方向；对于任一血管段中心线，具体步骤为：

6.1)找到其上的首尾两点P₁(x₁，y₁)、P₂(x₂，y₂)，并计算连接首尾两点所得的直线

6.2)根据点到直线的距离公式找到血管段中心线上到直线l距离最大的点P_max(x₃，y₃)，相应的距离为

其中A、B和C是将直线l的方程写为一般式Ax+By+C＝0时对应的系数；

6.3)如果距离D_max＜D_T，则将该直线作为该血管段中心线的近似，该血管段中心线处理完毕，其中D_T是预先给定的阈值；如果D_max≥D_T，则以P_max为分割点，将该血管中心线分为P₁P_max和P₂P_max两段，并分别对两段分别进行步骤6.1)～6.3)的处理；当所有分段都处理完毕时，依次连接各个分割点形成折线，作为该血管段中心线的近似。

Images