基于深度学习和二维高斯拟合的视网膜血管管径测量方法

刘海坤，王健，杨嵩，吴骏，尹昌顺，张凯，张震，裴新然，吴帅

1.天津工业大学电子与信息工程学院，天津300387；2.中国人民解放军空军93756部队教研部电子教研室，天津300131

前言

人体的眼底血管和心脑血管在生理解剖学上有共同的结构特性［1］，并且眼底血管是可以直观观察的血管。大量的临床经验表明，视网膜图像中血管管径的异常变化可以作为心脑血管疾病的诊断依据。例如，高血压视网膜病变的早期特征包括眼底动脉血管管径变窄［2］，中风和心脏病的早期特征则包括视网膜静脉变窄［3］。以上视网膜图像中血管管径的异常变化都出现在疾病初期，越早地发现血管的异常变化越能及早地确诊并治疗疾病，增加痊愈的概率。因此，对眼底图像血管参数进行测量在临床上具有重要的意义。

随着图像处理方法的发展，已有不少的科研人员对眼底血管管径测量问题进行了研究，目前现有的血管管径测量方法大体上可以分为血管边界法和横截面拟合法两类。

血管边界法是通过图像增强和分割的方法提取出完整的血管网络，获得血管边界，然后通过形态学方法提取血管的中心线，在骨架线上沿血管方向的垂直方向扫描，扫描到血管边界停止，计算边界两像素点之间的距离即为血管管径值。例如王淑珍等［4］通过对血管图像的一系列预处理和边缘检测方法测量血管；陆建荣等［5］利用血管的特性使用局部拟合血管中心线，在获得血管方向，利用方向和血管的中心线得到与血管边缘的两个交点，最后计算两点间距离得到血管管径；姚畅等［6］采用与王淑珍和陆建荣相同的方法，基于血管骨架线的位置信息和方向信息，利用血管周围的像素采用改进的定向对比度方法，得到垂直于血管骨架线方向的边界坐标求得血管管径。

横截面拟合法是在通过分析血管的横截面像素灰度分布情况，根据分布函数来拟合血管横截面的像素灰度求得血管管径。例如Zhou等［7］通过分析血管横截面的灰度分布，设计了一维高斯函数来拟合血管，找到血管边界值求得血管管径；Li等［8］根据不同血管的横截面的灰度分布特性，设计了不同的高斯函数来分别拟合不同类型血管，以此提高了管径测量的准确率。

根据国内外研究基础总结出的现有血管管径测量的方法可以看出，精确地分割血管是得到视网膜血管管径的前提和基础，目前自动识别视网膜血管的方法大体分为有监督和无监督的分类方法。有监督的方法是在有标准图像的基础上提取特征训练分类器，无监督的方法是在无标准图像的基础上设定判别准则。监督方法的性能往往高于非监督方法，但是传统的监督方法也有一些缺点，比如需要手动设计特征提取血管的方法，特征的设计必须要在选取分类器前定义好，对于特征的确定是一件费时费力的事，它需要设计者对所研究的对象很了解并且经验丰富，另外还有一大问题是需要根据不用的数据库需要重新选择设计特征，没有通用性。

深度学习是近年来人工神经网络领域的一项重大突破，通过组合浅层特征来形成深层特征，达到发现数据分布式特征的目的。和传统的监督方法相比，深度学习让计算机从观测数据中学习，根据学习结果自行解决问题。Liskowski等［9］从大图中提取图像块进行数据扩充，利用深度神经网络进行视网膜血管分割。Fu等［10］提出基于卷积神经网络（Convolutional Neural Networks,CNN）和条件随机场（Conditional Random Field,CRF）的视网膜图像血管分割方法，该方法将血管分割作为边界检测问题处理，利用CNN产生分割概率图，再结合CRF得到二值分割结果。全卷积神经网络（Fully Convolutional Network,FCN）［11］作为深度学习的一支，是基于语义级别图像分割问题提出的。以Ground Truth作为监督信息训练网络，让网络作像素级别的预测，使图像级别的分类进一步延伸到像素级别的分类。U-net［12］模型是基于FCN的语义分割网络，适合医学图像分割。该网络采用编码器-解码器结构，利用编码器逐渐减少池化层的空间维度，利用解码器逐步修复物体的细节和空间维度。本文在U-net网络部分结构的基础上，通过减少池化层来减少网络的深度，使网络更适合于学习视网膜血管的特征，同时，为了使网络更好地学习目标的细节，还将网络的血管特征层和血管位置层相连接，以学习的方式自动获取每个特征通道的重要程度，并依照重要程度自适应的调整各通道的特征响应，方便后续网络的特征提取，提高网络模型分割结果的精度。

1 基于改进型FCN网络的血管分割

1.1 深度可分离卷积

卷积层的主要作用是特征提取，标准的卷积将三维的卷积核作用在一组特征图上，需要同时学习空间上的相关性和通道间的相关性。卷积层的主要作用是特征提取，标准的卷积将三维的卷积核作用在一组特征图上，需要同时学习空间上的相关性和通道间的相关性。深度可分离卷积（Depthwise Separable convolution）在执行空间卷积的同时，保持通道之间分离，然后按照深度方向进行卷积。如图1所示，首先执行深度卷积（Depthwise convolution），即输入的每个通道独立执行空间卷积，然后执行1×1卷积（Pointwise convolution），将深度卷积的通道输出映射到新的通道空间。

标准卷积同时学习空间信息及通道间的相关性，然后对输出进行非线性激活；深度可分离卷积首先进行深度卷积，增加了网络的宽度，使得特征提取更加丰富，然后直接进行1×1卷积，对1×1卷积结果进行非线性激活。

假设有一个3×3大小的卷积层，其输入通道为16、输出通道为32。具体为，32个3×3大小的卷积核会遍历16个通道中的每个数据，从而产生16×32=512个特征图谱。进而通过叠加每个输入通道对应的特征图谱后融合得到1个特征图谱。最后可得到所需的32个输出通道。针对这个例子应用深度可分离卷积，用1个3×3大小的卷积核遍历16通道的数据，得到了16个特征图谱。在融合操作之前，接着用32个1×1大小的卷积核遍历这16个特征图谱，进行相加融合。这个过程使用了16×3×3+16×32×1×1=656个参数，远少于上面的16×32×3×3=4 608个参数。这个例子就是深度可分离卷积的具体操作，其中上面的深度乘数（Depth multiplier）设为1，这也是目前这类网络层的通用参数。当通道数量更多时，减少的参数量会更大。因此，深度可分离卷积不但能够拓展网络宽度，而且在一定程度上减少了参数量。

图1 深度可分离卷积示意图Fig.1 Schematic diagram of depthwise separable convolution

1.2 改进型FCN网络

本文将深度可分离卷积应用于FCN网络，对于深度可分离卷积，首先执行深度卷积，增加网络宽度，然后进行1×1卷积，融合通道信息；采用merge层将下采样层后的血管特征图和上采样后血管位置标注图叠加，以学习的方式调整特征参数，以实现更好的分割血管。

网络结构如图2所示，网络的左半部分包含两组Separable卷积层，Depthwise卷积核大小均为3×3，在每层Separable卷积后使用ReLu函数进行激活，每组Separable卷积后连接步长为2的2×2最大值池化。为了减少特征信息的丢失，每次降采样后都将通道数量加倍。网络的右半部分包含两组反卷积层，核大小为3×3。由于网络浅层提取的特征包含许多细节，深层的特征更加抽象，对浅层特征通道进行加权，以学习的方式调整不同特征通道的权重，增强有用特征，抑制不重要的成分，然后再将上采样的结果与通道加权后的浅层高分辨率信息进行连接。每组连接后再进行两次Separable卷积，Separable卷积通道的数量相比上采样的通道数量进行了减半处理，Depthwise卷积核大小均为3×3，同样在每层Separable卷积后使用ReLu函数进行激活。最后一层，使用卷积核大小为1×1的标准卷积将32个feature map映射为2个feature map，实现眼底图像血管以及背景二分类。

图2 网络结构图Fig.2 Network structure diagram

1.3 血管分割

1.3.1 图像预处理 为了使设计的网络能更好地学习血管的特征，需要对图像进行一系列的预处理，实现对血管的突出和对噪声的抑制。首先使用归一化操作。之后采用对比度受限的自适应直方图均衡化（Contrast Limited Adaptive Histogram Equalization,CLAHE）对归一化的图像进行处理，CLAHE同时具有自适应直方图均衡化和对比度限幅的优点，能提高视网膜图像血管的对比度和清晰度。最后对CLAHE的结果进行Gamma校正，提高图像的动态范围，实现对比度拉伸，增强图像对比度。原始图像预处理结果如图3所示。

1.3.2 数据增强 通过对血管特性的分析可知，血管的尺度不一，若是将整幅眼底图输入进网络进行特征的提取，会忽略毛细血管的提取，所以将整幅眼底平均分成图像块。基于图像块进行血管特征的学习，能更充分地学习血管规律不会错分毛细血管，因此，本文采用了图像块对网络进行训练。以DRIVE库中20幅图像为样本，每幅图像被裁剪成9 500个分辨率为48×48的图像块，20幅彩色眼底图共有190 000幅图像块。图像的裁剪解决了医学图像数据量不足的问题，尽管可能图像块重叠，即不同的图像块可能包含原始图像的相同部分，但不执行进一步的数据增强。本方法随机的将数据集的90%用于训练（171 000个图像块），剩余的10%用于验证算法（19 000个图像块）。训练集的血管金标准图也按照以上方式进行裁剪，以便神经网络在训练后进行验证使用。对彩色眼底图提取绿色通道后裁剪图和专家分割的血管裁剪图如图4所示，从图4a和图4b可以看到眼底图的裁剪图包括了含血管和不含血管的部分以及其对应的血管分割图。

图3 原始图像预处理结果Fig.3 Preprocess results of the original image

图4 图像数据集Fig.4 Image dataset

1.3.3 分割结果 图5给出了对DRIVE库测试集中的20幅眼底图的分割对比结果，第1列为本文网络分割的血管结果图，第2列为DRIVE库测试集中第一位专家分割的血管结果图，第3列为DRIVE库测试集中第2位专家分割的血管结果图。本文的分割结果和专家二的分割结果都和专家一的分割结果进行对比。从对比图中可以看出，本文的分割结果能完整的分割出血管网络，包括一些细小血管都可以被分割。和专家二的分割结果对比可以看出，一些专家二没有分割出的微血管本文方法也能正确分割。和专家一的分割结果对比可以看出，本文方法只有很少的微血管没有分割出。从以上的分析可以看出，本文可以解决现存的微小血管不能分割的问题。由于DRIVE库中的一些图存在出血点和硬渗等病变，从分割结果可以看出本文方法对于存在病变的眼底血管图也可以很好地分割出血管，解决了不能完整分割存在病变血管的问题。

图5 分割结果对比图Fig.5 Comparison of segmentation results

2 视网膜血管的管径测量

2.1 视网膜血管的中心线的提取

若要对血管进行管径测量，就要求得血管的方向，人眼可以很好地判断血管走向，但是机器是无法识别血管的结构和走向的，所以要先求取血管的拓扑结构，再根据拓扑结构就能方便求得血管的走向。血管中心线就是血管的拓扑结构，提取血管中心线的方法包括数学形态学提取中心线法、细化算法和中轴变换的方法。由于中轴变换和细化算法需要大量的计算，例如中轴变换就需要计算血管边缘上的每个点到每个血管区域内部点的距离，所以使用数学形态学方法提取血管中心线。数学形态学细化方法提取血管骨架，结果如图6所示。

图6 中心线提取结果Fig.6 Results of centerline extraction

2.2 视网膜血管方向的获取

由于血管曲率较小，可以将血管假设为分段的线性直线段。本文利用最小二乘法将血管段拟合为直线，求出血管段的斜率。最小二乘法是将一组符合y=a+bx关系的测量数据，用计算的方法求出最佳的参数a和b。

在研究两个变量(x,y)的对应关系时，通常可以得到一系列成对的数据(x1,y1、x2,y2,...,xn,yn)，把这些数据描绘在x-y坐标系中，若发现这些点在一条直线附近，则可设直线方程的表达式（1）为:

要根据测量数据求出最佳的a和b，令：

D对a和b分别求一阶偏导数和二阶偏导数，满足最小值条件，令一阶偏导数为零，引入平均值后解得血管斜率为：

由上述最小二乘法的原理，根据血管段骨架线上像素点的坐标，就可以求出当前血管段的斜率b。设当前像素点为Pi，则当前像素点Pi的方向如式（4）为：

2.3 基于二维高斯拟合的血管方向及中心线矫正

由于视网膜血管可能会出现弯曲的状态，采用最小二乘法拟合直线仍然存在一定的误差，如图7所示，在血管弯曲处求得的方向和血管实际走向存在一定的误差。

图7 传统方法求取的血管方向Fig.7 Blood vessel direction obtained by traditional method

对管径测量来说，形状特性是非常重要的，血管中心线恰能准确的表征血管的形状特性，所以若要准确测量血管的管径值，就要获得准确的血管中心线。利用数学形态学方法提取血管中心线，结果如图8所示。从图中黑色方框可以看出，由于形态学骨架线提取算法对所采用的结构元素有很强的依赖性，骨架线提取结果随着结构元素的结构和大小的改变而改变，因此在血管分叉处或交叉处提取的中心线的形状特性偏离血管的形态特征。

图8 传统方法求取的血管骨架线Fig.8 Blood vessel skeleton lines obtained by traditional method

针对以上两个问题，根据血管横截面灰度呈高斯分布的特性，根据血管初始的方向和中心线使用二维高斯函数对血管拟合，得到校正后准确的血管方向和位置信息，最后在血管分割结果基础上测量血管管径。

常用的二维高斯函数表达式为［13］：

在公式（1）中不包含血管的方向信息，为了拟合出血管方向信息，本文提出包含血管方向信息的二维高斯函数，数学表达式如式（6）、（7）所示：

其中(x,y)为校正前的中心线坐标，(x0,y0)为校正后的中心线坐标，θ为校正后的血管方向，σ为高斯函数的均方差，A为高斯函数的峰值，对于高斯函数的未知参数求解过程如下：

首先，对公式（8）两侧取对数并乘f，可得：

然后，假设在血管骨架线图像中共含i个候选像素点，其中(xj,yj)为候选像素点中的第j个像素点坐标，f(xj,yj)为第j个点在图像中的灰度值。将所有候选点像素带入二维高斯模型中，并对a、b、c、d、e进行求解：

最后，根据求得a、b、c、d、e的值计算中心线坐标及血管方向：

（1）中心线坐标 (x0,y0)的求解，根据公式（10）推导可得：

（2）血管方向θ的求解，根据公式（11）推导可得：

对在形态学细化后的骨架线上点A(xi,yi)的方向及中心线位置的校正如图9所示，其中l为血管骨架线，A点为骨架线的原始点，利用Hessian矩阵的多尺度增强滤波器得到A点的最大输出响应方向，即血管初始走向，再将点A的方向及位置信息输入进二维高斯函数后，利用血管横截面上的灰度值呈高斯分布的特点，对点A进行方向及中心线位置的校正，得到校正后的点B(xi,yi)及其方向。血管横截面灰度值特征函数如图10所示，其中B为校正后的血管中心点，从图中可以看出校正后点B的灰度值为血管横截面的最大值。

图9 中心线及方向校正示意图Fig.9 Diagram of centerline and direction correction

图10 血管横截面灰度值特征函数Fig.10 Vessel cross-section grayscale value characteristic function

最后，利用上述2个高斯模型参数对血管进行高斯拟合，得到校正后的血管中心线及方向，它们与校正前血管中心线及方向对比图如图11所示。图11a为校正后血管中心线图，图11b为校正前的血管中心线图。图11a和图11b相比，可以看出分叉处的血管中心线走向有明显的改善，和血管实际走向符合且更平滑。图11c为校正后的垂直血管走向图，图11d为校正前的垂直血管走向图，黑色箭头代表实际的垂直血管走向。通过对垂直血管走向的观察可以看出校正后的血管方向更符合血管的实际走向，校正前的血管方向和血管实际走向存在一定偏差。此处使用血管方向的垂直方向作为图示是为了能更清晰地观察血管走向。由图11可以看出，校正后得到的新的血管中心线和方向更准确，更具实用性。

2.4 血管管径测量

在血管分割结果基础上测量血管管径。管径测量包括3步，首先利用形态学细化得到血管骨架线，然后使用二维高斯函数对骨架线上的初始点的方向信息和位置信息进行校正，得到准确的血管方向和中心线位置，最后沿垂直于血管的方向扫描找到血管边界点，计算两边界点之间的距离获得血管的像素级宽度。

图11 校正后的血管中心线及方向和校正前的血管中心线及方向对比Fig.11 Comparison of vessel centerline and direction before and after correction

在血管的垂直方向扫描以定位血管的边界点，中心线上各点与X轴的夹角是一个重要参数。血管管径测量图如图12所示，其中图12a为血管管径测量图，θ为中心线上一点Pi的血管方向，即可得到当前点与X轴的夹角为θ，图中虚线为血管的边缘。为了找到血管垂直方向上的血管边界点Q1、Q2，以Pi为起点，以1个像素为步长，沿着与血管垂直的方向分别向血管两侧扫描，找到血管两侧的边界点Q1、Q2。计算Q1、Q2两点的欧式距离即为血管骨架上Pi点处的管径，图12b为视网膜图像中局部部分血管管径测量示意图。

图12 血管管径测量图Fig.12 Diagrams of vessel diameter measurement

3 实验结果分析

3.1 血管分割性能分析

采用准确度（Accuracy）、血管分割错误率（False Positive Rate）以及血管分割正确率（True Positive Rate）3个参数对本文的血管分割方法性能进行评价［14］。ACC是识别彩色眼底图像的像素为血管或是背景的正确识别率；FPR是将背景或噪声像素识别为血管的假血管识别率；TPR则是将血管像素正确识别为血管的真血管识别率。其公式如下表示：

其中，TP表示本文方法将实际是血管的像素识别为血管，FP表示本文将实际不是血管的像素错误识别为血管，FN表示本文将实际是血管的像素错误识别为非血管，TN表示本文将实际不是血管的像素识别为非血管。

本文方法得到的2个参数和专家2的血管分割指标值进行对比，以专家1的分割图作为金标准，指标对比表如表1所示。

表1 本文分割方法参数和专家2分割参数对比Tab.1 Comparison of the segmentation results of the proposed method and expert 2

从表1可以看出，本文的FPR和ACC两个参数均好于专家二分割的参数，TPR仅仅略低于专家二分割的参数。说明了本文的分割方法的性能已经达到了医生的分割水平。

为了能更好地说明本文方法的性能，将本文方法的参数和近年来代表性眼底血管分割方法的参数进行对比，指标对比表如表2所示。

表2 本文分割方法参数和近年分割方法参数对比Tab.2 Comparison of the segmentation results of the proposed method and other segmentation methods in recent years

从表2数据可以看出，本文方法的FPR和ACC两个参数均优于其他方法的参数，而TPR仅仅比文献［15］方法和文献［9］方法的参数值略低了一些。从总体来看，本文方法的分割性能不仅高于专家的分割性能，还能高于其它的血管分割方法，代表了本文所提出的方法具有高水平的分割性能。

3.2 血管管径测量性能分析

由于REVIEW数据库中有公开的眼底图像血管片段以及对应的管径测量数据，为了对比说明，对此数据库提供的血管片段进行管径测量。分别对REVIEW数据库中的3个图像集进行血管管径平均值和管径标准偏差平均值计算，并分别与不同方法测量结果进行对比，对比结果如表3、表4、表5所示。

从表3中可以看出，本文方法测量的管径平均值与金标准数据的差值大小为0.1，且能够达到最低的标准偏差0.11，从而验证了本文方法管径测量的稳定性。从表4、表5可以看出，由于VDIS图集中存在中央反射的眼底图像，CLRIS中存在病变的眼底图像，通过对VDIS和CLRIS数据库进行管径测量并与不同方法对比，本文方法测量的血管管径与血管管径金标准很接近，表明本方法能够有效测量存在病变的眼底图像血管管径。综合以上分析，验证了本文方法的准确性和稳定性。

4 结论

本文提出一种基于深度学习和二维高斯拟合的视网膜血管管径自动测量方法。首先将深度可分离卷积和通道特征叠加引用于全卷积神经网络中对血管网络进行分割。然后通过分析传统方法所提取的血管中心线和方向与实际的血管中心线和方向存在偏差，所以根据血管横截面灰度值呈高斯分布的特点利用改进的二维高斯函数拟合血管，校正血管的中心线和方向，最后进行血管管径和动静脉管径比值的测量。对REVIEW库中的KPIS图集、VDIS图集进行血管管径测量，通过和其他代表性方法的测量数据对比，本文的血管管径均值和金标准差距最小且标准偏差最小，证明了本文方法测量的结果准确性和稳定性。

表3 不同方法对KPIS图集管径测量数据对比（像素）Tab.3 Comparison of different methods for vessel diameter measurement in KPIS set(pixel)

表5 不同方法对VDIS图集管径测量数据对比（像素）Tab.5 Comparison of different methods for vessel diameter measurement in VDIS set(pixel)