糖尿病视网膜病变(Diabetic Retinopathy, DR)作为糖尿病的主要并发症之一, 其发病率逐年增加, 严重危害到人们的视觉质量, 是糖尿病患者致盲的主要原因. 临床上将已确诊为糖尿病, 且尚未出现明显眼底病变之前的时期称为DR临床前期. 而DR一旦由临床前期进展至临床期, 留给医生的时间和手段就会很有限, 而且预后效果总体不好. 如果医生能在DR临床前期进展至临床期的易发阶段进行有效的诊断和治疗, 对于改善患者视功能预后、提高患者生活质量具有十分重要的意义^[1]. 据调查显示, 糖尿病病程在10年以内DR患病率为7%, 15年以上者为63%, 30年以上为95%^[2]. 据此可知, 糖尿病病程在10年以上的DR临床前期患者更容易病变到临床期. 因此选取糖尿病病程在10年以上, 且尚未出现明显眼底病变的患者作为研究对象更具有参考意义.

眼底自发荧光影像作为近二十年来开展的一项技术, 是一种无创的检查手段, 广泛应用于视网膜疾病的诊断. 目前, 在医学图像模式识别领域, 对于眼底荧光图像的相关研究还很少, 而针对DR临床前期这一阶段的研究更是几乎没有. DR临床前期和正常的眼底自发荧光图像如图1所示, 由于DR临床前期的患者并没有出现明显眼底病变, 所以在视觉感官上, 临床前期和正常的眼底图像并没有明显的差异. 如果仅凭医生肉眼观察, 几乎是无法识别的. 但是在没有检眼可视的眼底病变的DR临床前期, 极有可能已经存在着神经和微血管两个方面的病变. 对此, 可大胆假设, 两类眼底自发荧光图像存在着由神经和微血管早期病变引起, 但视觉效果并不明显的纹理差异. 因此本研究尝试通过当前计算机视觉领域的一些纹理特征算法去建立DR临床前期的眼底自发荧光图像的识别模型.

纹理是一种重要的视觉线索, 是图像中普遍存在但又难以描述的特征^[3]. 纹理分析是计算机视觉和模式识别领域不可或缺的一部分. 局部二值模式(Local Binary Pattern, LBP)和灰度共生矩阵(Gray-Level Co-occurrence Matric, GLCM)是目前图像模式识别领域较为常见的两种算法, 均属于纹理分析方法中的统计分析方法. GLCM算法由于其不受分析对象的制约, 能很好地反映图像的空间灰度情况, 表现图像的纹理特征, 所以得到广泛应用. LBP在纹理分类上效果显著, 一经提出, 便在医学图像^[4]和人脸识别邻域^[5]取得了不错的应用.

1 实验数据和方法 1.1 数据来源

本研究的所有样本均于2018年2月至5月在瑞安市人民医院眼科门诊所收集, 其中受检者均进行前置镜、裸眼视力或矫正视力检查, 而且受检者需满足以下条件: (1)视力或矫正视力≥0.8且屈光度≤±6.0 D; 曾行白内障摘除人工晶状体植入手术者, 屈光度≤一2.0 D且手术前无高度近视病史; (2)裂隙灯显微镜检查屈光间质透明, 50岁以上者允许晶状体密度增加; 曾行白内障摘除人工晶状体植入手术者, 人工晶状体位置正常, 后囊无明显混浊; (3) 取得所有受检者的知情同意后, 检查双眼眼压后, 用0.5%复方托吡卡胺滴眼液散瞳, 瞳孔散大至7～8 mm, 由同一临床经验丰富的眼底病专科医师行前置镜检查, 排除视神经或视网膜疾病; (4)眼压≤21 mm Hg; (5)无青光眼家族史、高血压等全身性疾病; (6)排除曾接受屈光性手术者.

研究组: 确诊为2型糖尿病的患者, 共计56例, 112只眼睛(56只左眼, 56只右眼), 其中男性32例, 女性24例, 为DR临床前期组; 对照组: 为医院健康体检者, 经检查未患有糖尿病者及其他眼底病变, 共计39例(38只右眼, 35只左眼), 73只眼睛. 其中男性17例, 女性22例, 为正常组. 所有眼底自发荧光图像的原图大小为496×596像素, 在去除了无关的背景信息后, 大小均为496×496像素.

1.2 实验方法

医学图像模式识别过程^[6]一般如图2所示. 将所有数据按比例随机划分为测试集和训练集, 随后保持不变. 分类模型通过训练集数据建立, 测试集用来验证模型的有效性. 当只考虑验证特征的有效区分度时, 也可采用全部数据交叉验证的手段. 本研究中首先采用支持向量机(SVM)对全部数据进行十折交叉验证以便挑选出有效的特征提取算法及特征参数, 随后采取留出法将所有数据一次性划分训练集和测试集做进一步检验.

1.2.1 LBP特征提取

LBP是一种用来描述图片局部纹理特征的算子, 其核心思想就是以中心像素的灰度值作为阈值, 对其邻域灰度值进行阈值操作, 对于灰度值大于该阈值的像素点置1, 否则置0, 从而得到一组二进制数^[7]. LBP编码公式为:

$ T_{{\rm {LBP}}(X_c, Y_c)} = \sum\limits_{P = 0}^{P - 1} {s({g_p} - } {g_c}){2^p} $

(1)

$ s(x) = \left\{ {\begin{array}{*{20}{c}} 1 & {\left( {x \ge {\text{0}}} \right)} \\ 0 & {\left( {x < 0} \right)} \end{array} } \right.$

(2)

式中, g_c为中心像素(X_C, Y_C)的灰度值, g_p为邻域点P的灰度值. 对于灰度图像的任一空间给定位置(X_C, Y_C), 将其与其8邻域灰度值比较后得到一个二值集合, 将其有序排列后转化为一个8位无符号的二进制数. 半径为1, 邻域为8的圆形LBP算子计算示例如图3所示.

由上可知, 对于邻域为P的LBP算子存在着2^P种二进制模式, 而过多的二进制模式, 会使产生的直方图显的过于稀疏, 这样对纹理的提取、识别还是分类都是不利的. 为此, Ojala等人提出了一种叫“等价模式”^[8]的方法来对LBP模式进行降维. 据他们发现, 有一类模式在图像中出现频率极其高, 能够刻画90%以上的纹理特征, 这类模式就是等价模式, 而它们都有一个特性, 就是0->1或1->0 的跳变总次数小于或等于二, 如00000000, 00110000, 10000001等. 而判断一个二进制模式是否为等价模式, 最简单的办法就是将LBP值与其循环移动一位后的值进行按位相与, 计算得到的二进制数中1的个数, 若个数小于或等于2, 即为等价模式. 而除了等价模式以外的模式都归为一类, 称为非等价模式类, 如10010011(4次跳变), 01000001(3次跳变)等. 我们可以推出, 对于P邻域的LBP算子, 等价模式的个数为P×(P–1)+2个, 加上非等价模式类, 总计P×(P–1)+3种二进制模式.

上述的LBP特征具备灰度不变性, 但还不具备旋转不变的性质, 在整张图像转动的情况下, LBP值会发生变化. 为此, Pietikinen等人对LBP算子进行了扩展, 提出了具有旋转不变性的LBP算子^[8], 即不断旋转圆形邻域得到一系列初始定义的LBP值, 取其最小值作为该邻域的LBP值.P=8, R=1的LBP算子的旋转不变过程如图4所示, 通过旋转变化不但使LBP算子具备了旋转不变性, 同时也一定程度上起到了降维作用. 例如, 当P=8, R=1时, LBP值会有256种可能, 通过等价模式降维后会降低至59, 再通过旋转不变变化后, 维度会进一步降低至10.

1.2.2 GLCM特征提取

GLCM^[9]是一种通过研究灰度的空间相关特性来描述纹理的常用方法, 它通过计算图像中一定距离和一定方向上的两点灰度之间的相关性, 来反映图像在方向、间隔、变化幅度及快慢上的综合信息. 其定义为: 设一幅灰度图像大小为M×N, f(x, y)表示其在图像点(x, y)处的像素点所对应的灰度值. (x₁, y₁), (x₂, y₂)为图像上相邻距离为d, 角度为θ的点对, 其像素值分别为f(x₁, y₁)=i, f(x₂, y₂)=j, 而p(i, j; d, θ)则表示类似的像素点对在图像中出现的概率, 表示如下:

$ \begin{gathered} p(i, j;d, \theta ) = \# \{ ({x_1}, {y_1}) \in (M, N), \hfill \\ \;\;\;\;\;\;\;\;\;\;\;({x_2}, {y_2}) \in (M,N){\text{|}}d, \theta , f({x_1}, {y_1}) = i, f({x_2}, {y_2}) = j\} \hfill \\ \end{gathered} $

(3)

式中, #表示集合中的数目. 通过改变d和θ, 可以统计不同空间位置的像素对.

GLCM包含了图像的纹理信息, 但其包含的数据过多, 不适合用作纹理特征. 为此, 研究者通常在其基础上进一步获取二次统计量来作为图像的纹理特征. 而本文仅针对能量(Energy), 相关性(Correlation), 对比度(Contrast)和一致性(Homogeneity)4个特征进行讨论, 其计算公式如下所示:

$ Energy = \sum\limits_{i, j} {p{{(i, j)}^2}} $

(4)

$ Correlation = \sum\limits_{i, j} {\frac{{(i - \mu i)(j - \mu j)p(i, j)}}{{{\sigma _i}{\sigma _j}}}} $

(5)

$ Contrast = \sum\limits_{i, j} {|i - j{|^2}p(i, j)} $

(6)

$ Homogeneity = \sum\limits_{i, j} {\frac{{p(i, j)}}{{1 + |i - j|}}} $

(7)

1.2.3 支持向量机与K-折交叉检验

支持向量机^[10]是一种用于数据分类、模式识别以及回归分析的监督学习模型. 其主要思想就是寻求一个分类间隔最大的最优超平面对两类样本进行分类^[11], 尤其适用于小样本数据的分类. 本研究中采用的的是台湾大学林智仁教授开发的LIBSVM工具箱, 核函数选用径向基函数(Radial Basis Function, RBF). 而利用SVM建模的关键在于搜寻合适的惩罚因子c和RBF宽度参数gamma, 一般通过K-折交叉检验确定. 惩罚因子c表示对误差的容忍度, c过大会导致过拟合, 而过小则会导致欠拟合, 在模型准确率相等或相近的时候一般优先选择c小的那一组. gamma隐含地决定了数据映射到新的特征空间后的分布, gamma越大, 支持向量越少, 反之, 支持向量越多, 影响训练和预测的速度.

K-折交叉检验^[12]的做法是将实验样本均匀、随机地分K个子集, 每次取其中的K–1个子集作为训练集建立模型, 再将剩余的1个子集作为测试集用来验证模型的识别准确率. 不断更换训练集和测试集, 共计K次, 最后将这K次过程中的准确率作算术平均.

2 实验结果

本研究实验环境为Windows10操作系统, CPU2.6 GHZ, RAM8.00 GB, Matlab的版本为R2017a.

2.1 十折交叉检验结果

在实验中, 为了寻求对两组图片识别效果最好的特征集, 首先提取了185张眼底图片的LBP及GLCM特征. 再将这些特征单独或组合后放入SVM中进行十折交叉检验, 并通过十折交叉检验得到的准确率, 来选择最好的特征或特征组合.表1所示为实验中识别效果最好的几组特征及特征组合. 实验中提取的特征分别为: ① P=8, R=1, 通过等价模式降维且未经过旋转变化的59维LBP特征; ② P=8, R=1, 通过等价模式降维并经过旋转变化的10维LBP特征; ③ θ为0°, 45°, 90°, 135°, d=1的灰度共生矩阵的能量; ④ θ为0°, 45°, 90°, 135°, d=1的灰度共生矩阵的相关性; ⑤θ为0°, 45°, 90°, 135°, d=1的灰度共生矩阵的对比度; ⑥θ为0°, 45°, 90°, 135°, d=1的灰度共生矩阵的一致性.

由表1可知, LBP和GLCM对两组图像都有较高的识别率, 为此可认为两组图片在纹理上确实存在着肉眼无法识别的差异, 但相对GLCM而言, LBP算法的识别率更高, 效果更好, 无论是10-LBP还是59-LBP的准确率都在90%左右. 而且从特征的组合结果中可以看出, 两种算法中准确率最高的特征相互组合不但不会导致准确率上升, 反而会其下降. 对此可认为单独的LBP算法是这两种算法中识别这两组图像效果更好的算法.

2.2 留出法及R0C曲线

为了进一步验证LBP算法的识别效果, 将185张眼底荧光图像中两个类别以1: 1的比例划分到训练集和测试集中. 训练集中正样本37, 负样本56; 测试集中正样本36, 负样本56. 将训练集数据通过SVM建立的模型对测试集进行分类测试, 其中10-LBP算子的准确率为85.87%, AUC值为0.930, 59-LBP算子的准确率为88.12%, AUC值为0.943. 通过ROC曲线和AUC的值可知, 由LBP和SVM建立的模型具有较高的诊断价值. ROC曲线如图5所示.

从上述多组实验结果中可知, LBP算法和SVM建立的预测模型能有效的识别出处于DR易发阶段的糖尿病患者. 而这也证明, 处于DR临床前期的患者的眼底荧光图像与正常人的图像确实存在肉眼无法识别的纹理差异.

3 讨论

DR临床前期进展至临床期的病变高危期是一个相对模糊的时期, 而目前医学领域有不少人在研究DR临床前期的眼部变化时选择糖尿病病程10年以上者^[13,14]. 如果能及时、有效地识别和确定这一时期的患者, 对帮助降低DR患者的患病率及改善糖尿病患者的整体视觉质量都具有十分积极的意义. 而本文中提到的LBP算法对于这一时期的眼底自发荧光图像具有很好的识别效果. 相信通过数据源的不断补充和算法的改进, DR临床前期眼底自发荧光图像的辅助诊断技术在不久将来将应用于临床诊断. 而这对于糖尿病患者和相关医务人员来说, 无疑是一个巨大的福音.