基于图论的视网膜自动分层方法

唐宁，樊金宇，邢利娜，王晶，蒋天亮，李云耀，史国华

(中国科学院苏州生物医学工程技术研究所，苏州 215163)

引言

数十年来，眼底疾病和视力损伤的发生率正以惊人的速度上升。根据世界卫生组织统计，截至2021年，全球罹患不可恢复的眼底疾病或眼盲的患者超过22亿，相较2020年同期增长百分之三十。我国作为视力疾病的高发地区，患病或致盲的人数约占全球的五分之一。视觉障碍人群中，约半数患者的视力损伤是由未及时治疗所致，因此，对眼底疾病的早期预防和治疗至关重要。

对视网膜组织进行精确成像是眼底疾病诊断和治疗的必要前提。视网膜是位于眼球壁最内层的一层透明薄膜，最厚处仅为0.5 mm，从内至外共具有十层组织结构，依次为神经纤维层、神经节细胞层、内从状层、内核层、外从状层、外核层、内/外界膜、视网膜色素上皮层和脉络膜[1]。自1823年DeBuc等[2]将检眼镜用于直接观察视网膜开始，眼底相机、扫描激光检眼镜、血管荧光造影、光学相干层析等眼底成像方式相继出现，不同成像方式所获取的眼底图像各具优势，其中具备深度信息的光学相干层析技术对视网膜内部微小病变的观测最具优势。

光学相干层析技术(optical coherence tomography，OCT)自1991年被提出以来[3]，凭借其高分辨率、高灵敏度、非接触式断层成像等特点，在生物医学成像方面获得了广泛应用。随着研究的不断深入，OCT在眼科成像领域逐渐崭露头角，尤其在眼底疾病诊断中，光学相干层析技术已经逐渐成为不可或缺的一环。OCT的核心是采集不同深度组织背向散射光的干涉信号，对其进行解析，以获取组织内部微观结构信息。对于视网膜这种具有多层不均匀组织结构的成像对象，OCT可以根据各层组织散射光不同的干涉强度将其微观结构与微小病变清晰地显示出来，使得医生对眼底疾病的评估与诊断更为精确。

视网膜各层结构的厚度和形状是检测眼底疾病的重要指标。在早期使用OCT进行眼底疾病治疗的过程中，医生需要对视网膜图像进行手动分层及主观判读。由于OCT图像存在伪影噪声，且各层结构形状不规则、边界模糊、对比度低，人工分割的方式精度差、效率低。因此，高精度的自动视网膜分层十分重要。目前常见的视网膜分层方法有基于概率统计、活动轮廓、模式识别、动态规划以及图论的五种分层。基于概率统计的分层方法最早由Hee等[4]于1995年提出，通过层间像素灰度值的变化，对神经纤维层和视网膜层进行分割，后续经过Raja、Wang、Dodo、Ometto等[5-8]的改进与完善，实现了对OCT视网膜图像的九层自动分层，误差约为5个像素。基于活动轮廓的分层方法最早由Fernandez等[9]于2004年提出，用于视网膜黄斑区的病变研究，后续经过Kiaee、Amini、Chandrappoc等[10-12]的改进与完善，实现了对OCT视网膜图像的四层自动分层，误差约为3个像素。基于模式识别的分层方法最早由Fuller等[13]于2007年提出，利用支持向量机对3D视网膜OCT图像进行分割，但边界检测效果较差，后续经过Esmaeili、Kepp、Wang等[14-16]的改进与完善，最高可实现对OCT视网膜图像的八层自动分层，误差约为2个像素。基于动态规划的分层方法最早由Chiu等[17]于2010年提出，通过动态规划技术对视网膜8条边界进行提取，后续经过Maloca、Devalla、徐肃仲等[18-20]的改进与完善，实现了对OCT视网膜图像的八层自动分层，误差约为5个像素。基于图论的分层方法最早由Breger等[21]于2006年提出，通过图割算法与图最小算法对视网膜进行分割，后续经过Raja、Stankiewicz、Berenguer-Vidal、Liu等[5，22-24]的改进与完善，实现了对OCT视网膜图像的十层全自动分层，误差仅为2个像素。和其余四种分层方法相比，基于图论的分层方法所获取的分割信息最为全面，且误差小、精度高，能基本满足实时分层的需求。

本研究提出了一种基于图论的视网膜OCT图像自动分割方法，在将OCT图像表示为节点图后，根据像素强度对各边进行权值分配，再利用经典的Dijkstra算法[25]在限定区域内搜索最小加权路径，即可对视网膜的多层结构进行准确分割，成功实现对OCT视网膜图像的七层自动分层。

1 基本原理

本研究自动分层算法可主要分为如下三个步骤：首先，进行OCT图像权值分配；其次，进行视网膜端点的自动初始化；最后，利用Dijkstra最短路径算法完成视网膜的分割。自动分层算法的流程见图1。

图1 自动分层算法流程Fig.1 Automatic segmentation algorithm flow

1.1 OCT图像权值分配

为缩小搜索范围并避免视网膜分割时的跨层错误，采集视网膜OCT图像后，首先对其进行平直化操作。RPE层是视网膜信号最强的组织层之一，在OCT图像中表现为像素亮度高于其他组织。依托此先验知识，本研究设置A-scan中梯度最大的像素点为RPE层。然后找到RPE层的最低点，以此点为参照竖直移动其他视网膜信号，使视网膜被展平。经平直化操作后，对视网膜OCT图像进行权值计算。将OCT图像表示为节点图，每个节点对应一个像素，节点的连接为边，相互连接的边为路径，总权重最小的路径即位视网膜的层间分割线。权重分配示例见图2，三个连接节点之中13边的权重低于12边的权重，说明节点1偏好于节点3而非节点2，故13边被选为分割结果。只有为各边分配合适的权重，才能在后续确立最低加权路径时获得精确的分割结果。

图2 权值分配示意图(a).OCT图像节点示意图；(b).节点权值分配Fig.2 Weight assignment diagram(a).nodes of OCT images; (b).weight assignment of node

由于OCT图像中视网膜为水平结构，可利用OCT图像垂直方向上像素强度对视网膜各层进行区分，即根据梯度变化强度计算权值。权值计算方法见式(1)。

wab=2-(ga+gb)+Wmin

(1)

式中，wab是分配给连接节点a和b的边的权值，ga是节点a处图像的垂直梯度，gb是节点b处图像的垂直梯度，Wmin是设置的垂直边缘权值。式中将低权值指定给具有较大垂直梯度的节点对。在算法中，ga和gb被设置为0到1之间的值，Wmin=1×10-5。考虑到梯度的方向性，这些权值需要不断调整。例如，为了分割玻璃体NFL层边界，已知边界在较亮层上方呈现较暗层。相比之下，NFL-GCL层边界有一个从亮到暗的层变化。因此，在计算梯度提取适当层时，可以使用诸如[1，-1]和[-1，1]的图像边界。

1.2 视网膜端点的自动初始化

在为OCT图像各路径分配合适的权值后，需要选择路径的起始点和终止点进行后续的最低加权路径筛选。由于眼底视网膜图像是由多个分层结构所组成，如果手动选择分割层或估计相应层开始和结束的节点，效率低且人工选择无法选择到精准的边界位置。因此，本研究提出了一种自动初始化选择端点的方法。

该端点初始化算法基于眼底视网膜图像的分层条件，即分割的层需要延伸至OCT图像的全宽。后续所使用的Dijkstra算法目的是寻求最小加权的路径，因此，可以在图像的两侧添加一列额外的节点。这些节点具有任意强度值，在垂直方向上为它们组成的边指定最小权值Wmin。Wmin小于原始图像的邻接矩阵中任何非零权值，不仅可以保持新添加列中节点的连接性，而且初始端点可以在这些列的垂直方向上高速遍历。此外，新添加的列中可以任意指定开始和结束节点，初始端点能够在获取最小加权路径前沿着这些列自由移动。当完成图像分割后，即可删除在两侧添加的列，从而在无端点初始化错误的情况下进行准确的切割。

使用自动端点初始化技术分割图像的示例见图3。将两个垂直列以任意像素值和最小垂直边权值添加到原图像的任意一侧。为了达到最完整的遍历，将左上角和右下角节点初始化为起始点和中止点。根据像素强度为边指定权值，像素强度越低所赋予的权值越小。红线表示使用Dijkstra算法生成的分割最小路径。最后删除两侧添加的列，获取分割界限，实现在指定任意端点的情况下对视网膜各层边界进行准确的切割。

图3 自动端点初始化实例Fig.3 Automatic endpoint initialization

1.3 利用最短路径法完成视网膜的分割

确定了视网膜OCT图像各路径的起始点和终止点后，利用Dijkstra算法寻找最小加权路径。本研究所提出的自动分层算法是一种迭代算法，即分割出第一层边界后，通过分割结果限定新的搜索空间，从而进行对下一层边界的分割。按照视网膜各层边界的突出顺序，依次对玻璃体-神经纤维层(玻璃体-NFL)、神经纤维-神经节细胞层(NFL-GCL)、神经节细胞-内网状层(GCL-IPL)、内网状-内核层(IPL-INL)、内核-外网状层(INL-OPL)、外网状-外核层(OPL-ONL)、视细胞的内、外段(IS/OS)、视网膜色素上皮-脉络膜层(RPE-Choroid)进行分割。上一层的分割结果并不会影响后续层分割的精度，各层之间的分割结果是相互独立的。

玻璃体-NFL层边界和IS/OS层边界像素对比度高，是视网膜图像中两个最突出的边界层。利用Dijkstra算法从图像左上角起点至右下角终点搜寻最低权值路径。获得的切割路径即为玻璃体-NFL边界或IS/OS边界。与IS/OS边界不同，玻璃体-NFL层边界是视网膜的最顶层，上方无高反射率结构的存在，因此，可以通过检测高反射率结构是否存在分辨玻璃体-NFL层边界和IS/OS层边界。对OCT图像进行高斯滤波和二值化后，计算切割路径上方区域亮像素的像素数目占比，若占比超过0.025，说明存在高反射率结构，该分割路径为IS/OS层边界；反之，不存在高反射率结构，该分割路径为玻璃体-NFL层边界。当获得的分割路径是IS/OS层边界时，新的搜索空间限定于分割路径上方20像素的区域，以便进行玻璃体-NFL层边界的分割；当获得的分割路径是玻璃体-NFL层边界时，新的搜索空间限定于分割路径下方40像素的区域，以便进行IS/OS层边界的分割。

与玻璃体-NFL层边界和IS/OS层边界不同，视网膜其他层间边界并不明显，需要对搜索区域进行准确限制。在对OCT图像进行双边滤波、直方图均衡化、canny边缘检测及开闭运算后，可以通过先验知识将图像大致分为三个搜索区域，最上面的搜索区域为NFL层和GCL层，中间的搜索区域为GCL层、INL层和OPL层，最下面的搜索区域为IS/OS层和RPE层。首先，在最上面的搜索区域中对NFL-GCL层进行分割。基于扫描方向与被测眼确定颞侧位置，根据NFL层不同厚度确定分割路径搜索区域。由于NFL层薄厚差距在10个像素以内，较薄层搜索区域限制为玻璃体-NFL层边界下方10个像素内，较厚层搜索区域限制为玻璃体-NFL层边界下方9个像素内，利用Dijkstra算法在搜索区域内搜寻所得到的最低权值路径即为NFL-GCL层边界。接下来在中间的搜索区域中对IPL-INL层、INL-OPL层和OPL-ONL层进行分割。利用Dijkstra算法对这三个边界进行初步分割后，利用中央凹进一步细分搜索区域，实现边界精确分割。将玻璃体-NFL层边界到IS/OS层边界最小厚度列所在位置的两侧扩展30个像素，该区域中平均层厚度小于10像素的列，判定为中央凹区域。定位中央凹后，合并中央凹附近的分割路径，在保持搜索区域下边界位置不变的前提下，重新搜索最低权值路径，获得精确的IPL-INL、INL-OPL和OPL-ONL层边界。最后，在最下面的搜索区域中对RPE-脉络膜层边界进行分割。基于前面获得的IS/OS层边界，利用Dijkstra算法对该搜索区域直接搜索最低权值路径，即可得到RPE-脉络膜层边界。

2 实验结果与分析

2.1 视网膜图像分层结果

基于Dijkstra最短路算法对视网膜OCT图像的分层结果见图4。从上到下可以明显观察到所划分的眼底视网膜的玻璃体-NFL、NFL-GCL、GCL-IPL、IPL-INL、INL-OPL、OPL-ONL、IS/OS、RPE-Choroid。

图4 视网膜OCT图像分层结果Fig.4 Retinal OCT images segmentation

为了对自动分层算法的准确度进行评估，将视网膜图像自动分层的结果与手动分层的结果进行对比。将自动分层图与手动分层图相叠加，手工分层图进行虚化处理，保留手工划分线段。图5(a)中彩色分割线由自动分层算法得到，点状分割线由手动分层得到。将图中红框部分局部放大，两种分层方式所获得的视网膜边界基本贴合，肉眼观察不到明显差距，见图5(b)。

图5 自动分层结果与手动分层结果对比图Fig.5 Comparison of automatic segmentation and manual segmentation

随后对分层误差进行定量分析。选取5组不同视网膜的样本，每组采集100张OCT图像，对500张图像分别进行手工分层和自动分层处理。将处理后的图像导入MATLAB，利用图像矩阵信息获取各层分界线间的像素差，得到手工分层与自动分层的各层间距最大偏差、平均偏差和标准差。两种分层方法的偏差结果见表1。

表1 自动分层与人工分层的偏差结果Table 1 Deviation of automatic segmentation and manual segmentation

由表1可知，相较于手工分层，自动分层的最大偏差分布于4.12～7.97 μm，平均偏差分布于2.72～5.78 μm，标准差约为1.19 μm。ILM、NFL-GCL、INL-OPL层边界分割偏差较大，这是由于该部分曲率变化较大，尽管如此自动分层偏差均小于2个像素，满足了眼底疾病诊断和治疗的需求。此外，上述误差分析只取两种分层方法层边界最大间隔处的差值，部分差值过大的原因是自动分层方法采用了平滑过渡曲线，而人工分层会将层边界的细微凸起如正常形变的部分勾勒出来，但其值均在可接受范围内。除了视网膜层的高曲率处(如黄斑中心两侧)会导致相对大的偏差外，其余的视网膜层边界线基本都能重合在一起。

人工分层的准确度相对较高，但速度较慢，对一层视网膜层的边界划分通常需要1 min左右。而自动分层软件对单张OCT视网膜眼底图像的响应和处理时间通常能保持在3 s以内，且可以批量处理。由此可见，在分层精度达到满足诊断的情况下，利用自动分层软件，效率将显著提高。

2.2 病变视网膜图像的分层

为了对自动分层算法的实用性进行评估，利用自动分层算法对病变视网膜图像进行分析处理，分析病变视网膜的曲率对分割结果影响。首先，对中心浆液性脉络膜视网膜病变的OCT图像进行了自动分层，病变的OCT图像见图6(a)。中心浆液性脉络膜视网膜病变是最常见的眼底疾病之一，是指由于视网膜色素上皮屏障功能受损，脉络膜出现液体渗漏，导致视网膜色素上皮下积聚液体而引发的视觉扭曲与视力丧失。分层结果见图6(b)，可以明显观察到OPL-ONL层边界存在上凸形变(蓝色分割线)，符合中心浆液性脉络膜视网膜病变的特征，即视网膜神经上皮局部浆液性脱离。由此可见，本研究所提出的自动分割算法可以准确识别高曲率变化的病变视网膜层边缘，医生通过分割结果可以快速地针对病变进行分析与诊断。

图6 中心浆液性脉络膜视网膜病变Fig.6 Central serous chorioretinopathy

对糖尿病视网膜病变的OCT图像进行自动分层，病变的OCT图像见图7(a)。分层结果见图7(b)，在微弱病变区域，分割层次明显，可以准确地分辨各层视网膜结构；在严重病变区域，分割效果较差，这是由于视网膜间的脱离孔隙严重影响了自动分层结果。针对该情况，可以将视网膜间脱离孔隙的上边缘作为边缘进行识别分割，并以平滑曲线进行孔隙间的连接过渡，获取更为精确的病变视网膜层次结构，辅助医生对病变区域进行定位与分析。

图7 糖尿病视网膜病变Fig.7 Diabetic retinopathy

3 结论

本研究提出了一种视网膜OCT图像的自动分层方法，基于经典的Dijkstra最短路算法，在三个限定区域内搜索最低权值路径以获取视网膜各层结构边界。该方法能够精准地分割出视网膜七层边界，分层时间为3 s左右，分层偏差小于两个像素。该方法对于存在病变的视网膜OCT图像也能够提供较好的分层结果，在辅助诊断、病灶定位等方面具备一定的优势。尽管如此，本研究提出的自动分层方法仍有传统图论方法的局限性，即在形变严重的病灶区域无法准确地划分视网膜各层结构，为了能够在曲率变化大的视网膜结构处实现精准分割，后续可尝试机器学习的方法，对不同形变程度的眼底视网膜进行分层。