2. 北京儿童医院 眼科, 北京 100045
2. Ophthalmology Department, Beijing Children Hospital, Beijing 100045, China
青光眼是一种以视神经乳头结构改变为特征的进展性视神经病变,世界卫生组织将其列为全球第二大致盲眼病[1]。青光眼的早期症状不明显,诊断困难,计算机辅助青光眼诊断是解决这一问题的有效途径。
本文所设计的计算机辅助诊断算法基于两项重要的检查报告,即眼底照和视神经纤维层断层扫描(OCT)。在青光眼诊断中,评估眼底照中视神经乳头结构特征是诊断的重要环节。视神经乳头是视神经进入眼睛后部位置的入口,在眼底照中,视神经乳头处于一个明亮的椭圆形区域。在青光眼诊断中,一般通过计算视神经乳头中的视杯(optic cup)和视盘(optic disc)的半径比率,并利用下上鼻颞(ISNT)准则来判断。正常盘沿多数符合ISNT准则,即各象限盘沿宽度由宽到窄的顺序依次是下方、上方、鼻侧、颞侧。以鼻侧的盘沿宽度作为参照比较下方和上方盘沿是否有变窄及丢失,来判断是否有青光眼视神经损害。OCT报告中附有描述颞象限、上象限、鼻象限、下象限的视神经纤维层厚度曲线图,专业医师诊断时,以眼底照片的视神经乳头特征图和OCT报告中的视神经纤维层的形态特征同时作为诊断依据。
从2000年左右开始,多种计算机辅助诊断青光眼的算法被相继提出[2-7],但这些算法大多基于一种检查结果,而没有将多种检查结果进行融合。近年来,也有学者设计了融合多种检查结果的计算机辅助诊断青光眼的算法。如Silva等[8]提出一种融合OCT和视野检查结果的自动算法,Gopinath等[9]将光谱域OCT中的视神经纤维层厚度和毛细血管密度两者融合来预测青光眼。然而以上文献并没有将医生诊断青光眼时最为依赖的眼底照和OCT图像进行融合。这一方面不利于获得更准确的预测结果,另一方面也不利于医生理解和接受算法给出的预测结果,并将其应用于临床诊断。
本文首次将眼底图像特征和视神经纤维层厚度两种指标融合,提出了一种基于Dempster-Shafer(DS)证据推论[10]的多视图集成算法。利用本文方法可将两种不同的结果集成并给出预测结果,为证据冲突情况下如何通过融合两种指标作出最终决策提供解决思路。
1 集成学习 1.1 基本原理分类问题是机器学习中的一类基本问题,而集成学习是将多个学习器通过一定的规则融合起来的一种算法,其优势在于,对于同一个分类问题,通过集成可得到多个弱学习器,将多个弱学习器的学习结果结合起来用于最终预测,可以得到比单个学习器更高的预测精度。
同一特征下样本集的集成方法主要有Bagging、Boosting等经典算法,其应用范围极其广泛[11];多视图的集成学习方法主要方法有DS证据推论[10]、随机森林等方法。对多分类器的结果进行整合时,主要利用简单投票法、贝叶斯投票法及线性整合等。本文主要基于集成学习方法研究青光眼诊断问题。
1.2 分类算法医疗疾病诊断是常见的监督分类问题,本文采用支持向量机和逻辑回归两种分类算法。支持向量机在统计学习VC理论和结构风险最小化基础上,根据有限的样本信息在模型的学习精度和准确度上寻求最佳折衷,以期获得最好的推广能力,在解决小样本、非线性及高维模式识别中表现出许多特有的优势。它引入核函数,将低维数据映射到高维空间,从而解决低维数据在原始空间线性不可分的问题[12]。逻辑回归是一种广义的线性回归算法,常应用于疾病自动诊断、机器故障判断、信用违约预测等方面,它将数据拟合到logistic函数,输出预测概率,再根据预测概率判断所属类别[13]。
2 基于集成学习的算法 2.1 数据预处理本文针对集成学习采集有关于青光眼诊断的两种不同的检查数据,即眼底照片和OCT检查报告。这两种检查均为图像数据,由医生根据专业知识对眼底照中的视神经乳头视杯和视盘进行分割,标识出视盘线和视杯线,如图 1所示。由于视盘和视杯是近似椭圆,将它们平均等分成360个扇形区域,在每个扇形区域内,杯盘的半径比率通过计算视杯与视盘的面积比的平方根得到,即可得到样本i对应的360维杯盘比向量ai∈R360×1。
从OCT检查报告中提取视神经纤维层厚度曲线,如图 2所示。将曲线分割为360等份,计算每一等份的厚度均值组成360维的向量
${\mathit{\boldsymbol{b}}_i} \in {\mathit{\boldsymbol{R}}^{360 \times 1}} $ |
式中,bi是样本i的视神经纤维层厚度向量。bi与ai相对应,即每个样本均包含眼底照中的360维杯盘比向量和OCT中的视神经纤维层厚度360维向量数据。
2.2 特征提取在青光眼诊断中,通常关注眼底图像的杯盘比的大小以及杯盘比形态是否符合ISNT准则。根据青光眼疾病特点,需要提取杯盘比的两维特征平均杯盘比和下上鼻颞评分。首先利用局部平均法去除视杯及视盘轮廓线中的噪声,并利用样条插值和曲线标准化得到反映盘沿结构的视盘视杯半径比曲线;再通过与基准线作内积运算得到下上鼻颞评分后,与平均杯盘比组成杯盘比的特征向量,来分别刻画平均杯盘比和盘沿形态符合ISNT准则的程度[4]。
正常眼的视神经纤维层厚度曲线具有明显的双峰特征(图 2),其纤维层厚度要大于青光眼患者,并且视神经纤维层厚度的形态特征符合ISNT准则。故本文以视神经纤维层厚度向量的均值和下上鼻颞评分作为视神经纤维层厚度的二维特征,借鉴文献[4]的思想提出一种从OCT检查报告中提取视神经纤维层厚度尺度和形态特征的新方法。首先基于OCT检查报告提取视神经纤维层厚度曲线变化的360维向量,然后基于杯盘比特征提取的方式从视神经纤维层厚度曲线中提取两维特征,分别用于刻画平均视神经纤维层厚度及其符合规则双峰结构的程度。
视神经纤维层厚度二维特征提取过程如图 3所示。
视神经纤维层厚度二维特征计算方法如下。
1) 选取l个具有典型双峰结构特征的正常眼视神经纤维层厚度曲线,计算这l个样本的原始视神经纤维层厚度向量的平均值向量c1,其每个分量计算式为
$ \mathit{\boldsymbol{c}}_1^i = \frac{{\sum\limits_{j = 1}^l {\mathit{\boldsymbol{g}}_j^i} }}{l},i = 1,2, \cdots ,360 $ | (1) |
式中,g为正常人的视神经纤维层厚度向量。
2) 计算c1的平均值c
$ \bar c = \frac{{\sum\limits_{i = 1}^{360} {\mathit{\boldsymbol{c}}_1^i} }}{{360}} $ | (2) |
3) 计算视神经纤维层厚度的标准化基准向量c2,其每个分量式为
$ \mathit{\boldsymbol{c}}_2^i = \frac{{\mathit{\boldsymbol{c}}_1^i - \bar c}}{{\frac{1}{{360}}\sum\limits_{i = 1}^{360} {\left| {\mathit{\boldsymbol{c}}_1^i - \bar c} \right|} }},i = 1,2, \cdots ,360 $ | (3) |
4) 对于第i个样本,计算其视神经纤维层厚度向量bi各个分量的平均值,作为视神经纤维层厚度曲线的第一维特征xi2, 1
$ \mathit{\boldsymbol{x}}_i^{2,1} = \frac{{\sum\limits_{j = 1}^{360} {\mathit{\boldsymbol{b}}_i^j} }}{{360}} $ | (4) |
5) 对于第i个样本,计算其视神经纤维层厚度的标准化向量,其每个分量计算如式(5)
$ \mathit{\boldsymbol{\widetilde b}}_i^j = \frac{{\mathit{\boldsymbol{\widetilde b}}_i^j - \mathit{\boldsymbol{x}}_i^{2,1}}}{{\frac{1}{{360}}\sum\limits_{j = 1}^{360} {\left| {\mathit{\boldsymbol{\widetilde b}}_i^j - \mathit{\boldsymbol{x}}_i^{2,1}} \right|} }},i = 1,2, \cdots ,360 $ | (5) |
6) 对于第i个样本, 计算其标准化向量
$ \mathit{\boldsymbol{x}}_i^{2,2} = \mathit{\boldsymbol{\tilde b}}_i^{\rm{T}}{\mathit{\boldsymbol{c}}_2} $ | (6) |
最终,第i个样本的360维视神经纤维层厚度向量bi被转化为二维特征向量(xi2, 1, xi2, 2)。在医学诊断中,青光眼的这两项指标通常都较小。
2.3 基于多视图的青光眼分类算法本文利用DS证据推论[10]的组合方法解决两个视图分类结果不一致的问题。基于DS证据组合规则,利用logistic函数分别输出单个样本在两个视图下的预测概率,并比较两个预测概率大小来判断最终的预测类别;同时以logistic函数作为概率函数,输出对应样本的后验概率值。logistic函数形式如下
$ h(y)=\frac{1}{1+\mathrm{e}^{-y}} $ | (7) |
式中,h(y)为输出概率,y为输入函数值。根据logistic函数,输出样本类别为1的概率是pi,类别为0的概率是1-pi。
通过支持向量机分类模型最终输出样本对应的模型预测类别和决策函数值。支持向量机的分类超平面为wTx+b=0,其中w=(w1, w2, …, wN)是分类超平面的参数。则样本i的决策函数为
$ y^{i}=\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x}_{i}+b $ | (8) |
式中,yi为样本i的决策函数值,b为偏置,xi为输入向量参数。
样本i与分类超平面的欧式距离定义为
$ d_{i}=\frac{\left|\boldsymbol{w}^{\mathrm{T}} \boldsymbol{x}_{i}+b\right|}{\sqrt{\|\boldsymbol{w}\|^{2}}} $ | (9) |
在支持向量机中,通常认为样本点离分类超平面越远,分类预测的概率越高。对于支持向量机分类预测的概率输出,通过logistic函数将支持向量机所输出的决策函数值映射至[0, 1]区间,作为分类预测的概率输出[14]。
对于样本i的眼底照数据xi1和OCT数据xi2,基于眼底照视图和OCT视图下的训练集合,利用分类算法分别得到预测函数f1(x)、f2(x)。基于logistic函数输出样本i预测类别为1的概率,记为pji
$ p_j^i = \frac{1}{{1 + {{\rm{e}}^{ - {f_j}\left( {\mathit{\boldsymbol{x}}_i^i} \right)}}}},j = 1,2 $ | (10) |
则样本i预测类别为0的概率分别是1-p1i和1-p2i。
基于DS证据推论的集成学习算法将杯盘比特征和视神经纤维层厚度特征分别看作两个证据,根据组合规则融合基于二者的分类器输出的概率,输出预测类别为1和0的融合概率分别如式(11)、(12)所示
$ m_{1}^{i}=\frac{p_{1}^{i} p_{2}^{i}}{p_{1}^{i} p_{2}^{i}+\left(1-p_{1}^{i}\right)\left(1-p_{2}^{i}\right)} $ | (11) |
$ m_{0}^{i}=\frac{\left(1-p_{1}^{i}\right)\left(1-p_{2}^{i}\right)}{\left(1-p_{1}^{i}\right)\left(1-p_{2}^{i}\right)+p_{1}^{i} p_{2}^{i}} $ | (12) |
式中,m1i和m0i分别为样本i预测类别为1、0的概率,且m1i+m0i=1。比较m1i和m0i的大小,若m1i>m0i,则最终预测类别为1,反之为0。
基于DS证据推论的集成学习算法的基本流程如图 4所示。
本文数据来源于北京同仁眼科医院,均取自15天内同时做过眼底图像、视神经纤维层断层扫描检查的正常人和青光眼患者检查结果;共105个样本,包括59个正常眼和46个青光眼。采用留一法进行交叉验证,训练方式为:每次选取1个样本作为测试集,剩余104个样本作为训练集,所有的分类预测共循环105次,每次预测1个样本。计算105次预测准确率的平均值作为本文预测的准确率,并计算相应的特异性和敏感性;循环预测105次,分别计算平均准确率、平均敏感性和平均特异性。
敏感性(或真阳性率, Psen)用于表示诊断方法对疾病的识别程度,敏感性越高,漏诊率越低,其计算公式为
$ P_{\mathrm{sen}}=\frac{T_{\mathrm{p}}}{T_{\mathrm{p}}+F_{\mathrm{n}}} \times 100 \% $ | (13) |
特异性(或真阴性率, Pspe)用于表示未患病的确诊率,特异性越高,确诊概率越高,其计算公式为
$ P_{\mathrm{spe}}=\frac{T_{\mathrm{n}}}{T_{\mathrm{n}}+F_{\mathrm{p}}} \times 100 \% $ | (14) |
式中,Tp、Tn、Fp、Fn分别为真阳性、真阴性、假阳性、假阴性的样本个数。
受试者工作特征(ROC)曲线将假阳性率定义为X轴,真阳性率定义为Y轴。area under cure(AUC)是一种评价分类模型性能优劣的指标,被定义为ROC曲线与坐标轴围成的面积,AUC值越大,分类器性能越好。本文选择AUC值作为评价指标之一。
3.2 特征提取基于眼底图像和OCT报告提取360维的向量,计算向量的均值以及下上鼻颞评分作为二维特征,其二维特征分布如图 5和图 6所示。可以看出,本文提取的特征向量在二维空间内均呈线性不可分的特性。
在不同分类器下,用单个特征杯盘比或视神经纤维层厚度以及本文提出的集成学习算法融合两种特征进行分类预测,得到分类预测效果分别如表 1、表 2所示。以支持向量机作为基础分类器,以单个特征进行分类预测的准确率分别为87.62%、92.38%,而集成学习算法的准确率达到了95.24%;以逻辑回归算法作为基础分类器,集成学习算法的最终预测准确率达到了94.29%;同时本文算法预测的敏感性和特异性也高于其他两种方法。以上结果说明集成学习算法在线性和非线性分类器下的预测准确率、确诊率及漏诊率与传统算法相比都得到了较大的改善,利用本文算法可以得到更好的分类效果。
表 3为利用文献[3-7]的方法在本文数据集上预测的结果。图 7是分别以支持向量机、逻辑回归作为基础分类器,以本文算法进行分类的ROC曲线图,经过多视图的集成后,模型的AUC值分别为0.992、0.985。结合表 3可以看出,本文集成算法的AUC值和准确率均高于其他算法,并且以支持向量机作为基本分类器的AUC值达到了最高0.992。原因是本文算法的泛化能力较好,对于弱分类器的分类预测性能有着较强的提升能力,而其他算法只基于杯盘比或视神经纤维层厚度的特征进行分类预测,导致分类性能较弱。
基于医学专家知识,提出集成学习算法,对杯盘比和视神经纤维层厚度原始数据降维后分别提取出新的二维特征来代表尺度和形状信息,然后将杯盘比与视神经纤维层厚度这两种特征融合起来联合诊断青光眼。分别采用线性分类器和非线性分类器预测分类,对两个视图上的预测进行融合后,预测准确率均有较大的提高;在真实数据集上的实验结果也验证了本文算法较优的预测准确率、敏感性和特异性。
[1] |
KLEIN B E K, KLEIN R, SPONSEL W E, et al. Prevalence of glaucoma-the beaver dam eye study[J]. Ophthalmology, 1992, 99(10): 1499-1504. DOI:10.1016/S0161-6420(92)31774-9 |
[2] |
SALAM A A, AKRAM M U, WAZIR K, et al. Autonomous glaucoma detection from fundus image using cup to disc ratio and hybrid features[C]//2015 IEEE International Symposium on Signal Processing and Information Technology (ISSPIT). Abu Dhabi, 2015: 370-374.
|
[3] |
DAS P, NIRMALA S R, MEDHI J P. Diagnosis of glaucoma using CDR and NRR area in retina images[J]. Network Modeling Analysis in Health Informatics & Bioinformatics, 2016, 5(1): 1-14. |
[4] |
XU Y L, JIA X, HU M, et al. Feature extraction from optic disc and cup boundary lines in fundus images based on ISNT rule for glaucoma diagnosis[J]. Journal of Medical Imaging & Health Informatics, 2015, 5(8): 1833-1838. |
[5] |
BIZIOS D, HEIJL A, HOUGAARD J L, et al. Machine learning classifiers for glaucoma diagnosis based on classification of retinal nerve fibre layer thickness parameters measured by stratus OCT[J]. Acta Ophthalmologica, 2010, 88(1): 44-52. DOI:10.1111/aos.2010.88.issue-1 |
[6] |
GROUP A I F G S. Combining information from three anatomic regions in the diagnosis of glaucoma with time-domain optical coherence tomography[J]. Journal of Glaucoma, 2008, 23(3): 129-135. |
[7] |
LARROSA J M, POLO V, FERRERAS A, et al. Neural network analysis of different segmentation strategies of nerve fiber layer assessment for glaucoma diagnosis[J]. Journal of Glaucoma, 2015, 24(9): 672-678. DOI:10.1097/IJG.0000000000000071 |
[8] |
SILVA F R, VIDOTTI V G, CREMASCO F, et al. Sensitivity and specificity of machine learning classifiers for glaucoma diagnosis using spectral domain OCT and standard automated perimetry[J]. Arquivos Brasileiros De Oftalmologia, 2013, 76(3): 170-174. DOI:10.1590/S0004-27492013000300008 |
[9] |
GOPINATH K, SIVASWAMY J, MANSOORI T. Automatic glaucoma assessment from angio-OCT images[C]//IEEE 13th International Symposium on Biomedical Imaging (ISBI). Prague, 2016: 193-196.
|
[10] |
ZADEH L A. A simple view of the dempster-shafer theory of evidence and its implication for the rule of combination[J]. Ai Magazine, 1986, 7(2): 85-90. |
[11] |
ZHANG C, MA Y. Ensemble machine learning methods and applications[M]. New York: Springer Science & Business Media, 2012: 1-34.
|
[12] |
CORTES C, VAPNIK V. Support-vector networks[J]. Machine Learning, 1995, 20(3): 273-297. |
[13] |
COX D R. The regression analysis of binary sequences[J]. Journal of the Royal Statistical Society, 1958, 20(2): 215-242. |
[14] |
PLATT J C. Probabilities for SV machines[M]//SMOLA A J, BARTLETT P, SCHOLKOPF B. Advances in large margin classifiers. Cambridge: The MIT Press, 1999: 61-74.
|