高亚明 刘兆邦 陈 斌 李 铭 黄来剑
1(扬州大学信息工程学院 江苏 扬州 225127)2(中国科学院苏州生物医学工程技术研究所 江苏 苏州 215163)3(温州市人民医院 浙江 温州 325699)
近年来,我国人民的生活质量逐步提高,尿路结石的发病率也随之递增,南方发病率高于北方,约为22%~45%,个别省市甚至高达50%以上,且25%的患者需要住院治疗。显而易见,尿路结石已成为重要的公共卫生问题[1]。国内大多数研究表明,尿路结石的形成与很多因素有关,如地理环境、气候、饮食习惯等。本文主要研究体内结石中的单纯性结石,其成分以草酸钙最为多见,占结石的80%以上,无水尿酸次之[2]。目前利用传统的结石成分分析方法对两者进行区分需要借助仪器,步骤复杂且成本较高,无法实现术前体内无创检测。因此,本文通过CT图像,结合机器学习算法对两者的结构进行分析鉴别[3]。
就目前而言,基于机器学习的辅助诊断算法已经从纯粹的理论发展到了临床试验,在诸多方面如乳腺、肺结节等诊断上证实了其可行性与适用性[4],但其在结石成分鉴别上的研究较少,比如Perrot等[5]利用放射和机器学习能准确鉴别肾结石和静脉增生,却无法有效地在结石间进行鉴别。因此,针对体内单纯性结石特征,本文将辅助诊断算法用于其成分分析中,辅助医生进行诊断。
预处理主要运用图像增强技术,传统的图像增强技术大多是基于空间域对图像进行处理,本文利用灰度调整和图像插值方法。
灰度调整利用灰度直方图得到图像中的像素亮度分布情况,再通过均衡化、规范化处理,调整图像的清晰度。
本文首先将CT图像的普通像素值更改为医学专用的CT值,再根据结石组织CT值范围截取[100,1 500]HU以内的内容进行灰度调整[6],其中HU(Hounsfiled Unit)值表示组织对X射线的吸收程度。调整效果如图1所示。
图1 灰度调整
图像间插值在每帧图像中增加插值图像以缩小距离,并将插值图像和原始图像共同组成三维数据以提高三维图像的质量[7]。
本文利用立方插值算法在图像序列间进行超分辨率重建,将CT图像的像素间距统一为0.7 mm×0.7 mm×0.7 mm,如图2所示。
图2 图像间插值
当前医学图像分割正从手动或半自动向全自动发展[8]。本文将临床专家人工分割的结果作为金标准,利用3D Slicer软件对病灶组织手动勾勒感兴趣区域(ROI)[9]并对其裁剪得到更精准的ROI窗口[10],分割过程如图3所示。
图3 图像分割
由于体外的草酸钙结石质硬,呈环形或桑葚形,表面粗糙有刺,而无水尿酸结石质硬且表面光滑,呈圆形或卵圆形[11],故本文针对体内结石CT图像,分别提取在二维和三维上的灰度、形状和纹理特征[12],并对其进行对比分析。
医学图像中的灰度直方图考虑形成图像像素的强度,定义如下:由f(x,y)表示一幅灰度图像,(x,y)处的值表示该位置像素,一幅图有M个像素,分布在0~L-1灰度级,灰度直方图则统计每级灰度像素数目得到的统计图。
(1)
式中:i表示灰度级;L表示灰度级种类数(L≤256);mi表示灰度级为i的像素个数;M表示图像总像素个数。
本文利用灰度直方图,得到的灰度特征如下:
1)最大灰度值:分布在CT图像中的最大灰度级。
2)均值:图像灰度平均值。
(2)
3)方差:图像灰度在数值上的分布情况。
(3)
4)熵:直方图分布的均匀性。
(4)
5)倾斜度:直方图分布的不对称程度。
(5)
6)峰度:图像的灰度分布在接近均值时的大致状态。
(6)
形状特征为结石分析提供重要的结构信息,分为轮廓特征和区域特征。轮廓特征主要针对结石的外边界,本文得到基于轮廓的特征如下:
1)紧密度:衡量一个形状的紧致程度。
(7)
式中:A为形状面积;i1和i2为二阶矩:
式中:image(i,j)为CT图像;(i,j)处的值为CT图像中该位置的像素值;I为CT图像中所有像素点的横坐标集合;J为CT图像中所有像素点的纵坐标集合。
2)最长径:边界上相距最远的两个点之间的距离。
3)离心率:用焦点间的距离除以长轴的长度。
区域特征则针对整个结石区域,本文得到基于区域的特征包括二维区域中的面积和三维区域中的体积。
纹理特征反映结石表面的粗糙度、光滑性[13],提取方法分为统计分析、模型分析、结构分析和频谱分析。
本文采用基于统计的方法,主要包括:
1)灰度共生矩阵(GLCM):统计图像中不同灰度值的像素对出现的情况。
2)灰度游程矩阵(GLRLM):统计图像中具有相同灰度的直线区域(灰度、方向、长度)出现的情况。
3)灰度区域尺寸矩阵(GLSZM):统计图像中具有相同灰度的区域(灰度、大小)出现的情况。
4)邻域灰度差分矩阵(NGTDM):统计图像中相邻区域的灰度差异情况。
5)灰度差分统计(GRAY_DIFF):反映目标场景在连续时间点图像相减所构成的图像特征。
其中灰度共生矩阵(GLCM)最为经典,应用最广,由它提取出如下特征:
1)角二阶矩(energy),反映图像区域的均匀性或平滑性,表示为:
(8)
2)对比度(contrast),表示图像的纹理清晰度,表示为:
(9)
3)相关系数(correlation),反映矩阵行与列的线性相关程度,值越大图像区域灰度分布越均匀,表示为:
(10)
式中:μx、μy、σx、σy分别定义为:
4)熵(entrop),度量图像内容的随机性,表示为:
(11)
5)逆差分矩(uniformity),定义为:
(12)
6)同质化(homogeneity),定义为:
(13)
7)和平均(sum average),定义为:
(14)
式中:k=i+j。
所提特征如表1所示。
表1 所提特征
特征提取伪代码如下:
功能:特征提取。
输入:路径字符串1,路径字符串2。
输出:2D灰度特征,2D形状特征,2D纹理特征,3D灰度特征,3D形状特征,3D纹理特征。
从磁盘读取dcm医学图像数据文件
从磁盘读取nrrd医学图像标签文件
forx=3 topatient_num
获取dcm文件与nrrd文件的数量img_num
fori=1 toimg_num
按照InstanceID对dcm文件进行排序,保证与nrrd文件顺序一致
end
fori=1 toimg_num
查看当前帧是否存在标签
ifm1>0 then
对当前帧图像预处理
fory=1 tonum
对当前帧图像的每个连通区域进行分割,得到ROI
对ROI进行裁剪得到最大ROI,即CT
将所有图像Im和标签mm分别保存为三维矩阵O_CT3V、MSK3V
end
end
end
对CT提取二维灰度、形状、纹理特征
将所提特征归类为二维全局特征和内部特征
对O_CT3V,MSK3V进行裁剪
计算维度dim
ifdim>1
对三维图像预处理
fory=1 tonum
对每个三维连通区域进行分割,得到ROI
遍历得到最大的ROI,并存储为CT3V
end
forz=1 tosize(CT3V,3)
过滤,剔除没有标签的帧图像
end
对CT3V提取三维灰度、形状、纹理特征
将所提特征归类为三维全局特征和内部特征
else
只有一帧图像,直接将三维当作二维处理,得到特征
end
end
特征选择是对高维特征进行筛选以降低特征维度的过程[14],由于最小冗余最大相关(mRMR)算法综合考虑了相关性和冗余性,故本文使用mRMR算法。
mRMR算法通过计算特征之间和类标签之间的互信息来选出冗余性最小和相关性最大的前N个特征,即从最小冗余和最大相关出发,给出一种基于互信息的评价准则。其中:最大相关指的是特征与类标签信息之间的相关度最大;最小冗余指的是特征之间的冗余度最小[15]。
此外,将mRMR算法与相关系数(Spearma)法、Relief算法、SVM_RFE算法进行对比分析,从而验证其在本文中的适用性。
本文采集的CT图像数据是有限的,而支持向量机(SVM)分类器本身能较好地解决小样本、高维度等问题,故本文选择SVM分类器,针对有限样本进行训练。
在核函数的选择上,基于高斯径向基核函数的高维映射能力强、样本适用性好、多用于特征数量小(15维)、样本数量正常(119例)的情况,故本文使用基于高斯径向基核函数的支持向量机(RBF_SVM)分类器。具体实现则采用第三方LIBSVM工具包,利用其内置参数进行调参。
同时,将其与基于Linear核函数的SVM(Linear_SVM)分类器,随机森林(RF)分类器和Adaboost分类器进行对比分析,从而验证RBF_SVM分类器在本文中的适用性。
1)实验所用数据:
从温州市人民医院采集结石CT图像DICOM数据集,其中:草酸钙患者59例;无水尿酸患者60例;层内像素分辨率为512×512,层距为5 mm。
2)实验所用机器环境:
操作系统:Windows 10;CPU:Intel(R)Core(TM)i5-7500 CPU @ 3.40 GHz(3 401 MHz);内存:16 GB;编程软件:MATLAB R2018a。
1)分层10次10折交叉验证法。分层10折交叉验证首先将数据集分成10份,每份中类别之间的比例与整个数据集中的比例相同,轮番将其中9份作为训练集,1份作为测试集进行实验,每轮得到一个精度,最后对10轮结果取平均作为最后的模型精度。模型精度估计需要进行多次10折交叉验证,本文进行10次10折交叉验证,对得到的10个精度(ACC)求均值。
2)调参。参数设定直接影响算法性能,本文利用交叉验证法,结合ROC曲线下的面积(AUC)对各分类器参数设置步长,进行调参,以选取最优参数,优化分类模型。
3)性能度量。混淆矩阵是一个误差矩阵,由它可得许多度量指标,如表2所示。
表2 混淆矩阵
(15)
(16)
(17)
(18)
(19)
(20)
ROC曲线的横轴为假正例率:
(21)
ROC曲线的纵轴为真正例率:
(22)
根据分类器的预测概率迭代更新阈值,得到多个混淆矩阵对应的坐标点,依据坐标点画出ROC曲线进而计算AUC值。
如图4所示,首先对二维图像进行灰度调整,对三维图像序列进行插值,完成超分辨率重建,统一分辨率;其次,利用已标好病灶的NRRD标签文件进行图像分割,并裁剪得到更精准的ROI窗口;然后,再对ROI分别提取二维和三维上的灰度、形状、纹理特征,并通过特征选择算法进行筛选;最后将筛选得到的特征放入分类器进行模型训练,通过ACC、AUC等指标对比分析分类效果。
图4 实验步骤
1)特征提取算法对比。如表3和表4所示,各领域信息互补,综合训练模型的鉴别性能,在统一使用mRMR算法与RBF_SVM分类器的前提下,纹理特征比灰度和形状特征鉴别能力更强,2D特征略优于3D特征,而提取所有特征时,ACC值和AUC值最高,即最能反映病灶组织空间分布的异质性。
表3 特征提取平均ACC比较 %
表4 特征提取平均AUC比较 %
2)特征选择算法对比。对各特征选择算法进行对比分析,由表5可知,在统一提取所有特征与使用RBF_SVM分类器的情况下,mRMR综合考虑了特征之间的冗余性和特征类标之间的相关性,故其选择效果最好,验证了mRMR在本文中的适用性。
表5 特征选择对比 %
3)分类器对比。由表6可知,在统一提取所有特征与使用mRMR的情况下,RBF_SVM因其在小样本上的适用性较强而分类效果最好,验证RBF_SVM在本文中的适用性。
表6 分类器对比 %
4)最终辅助诊断算法评估。对结石CT图像预处理及分割后提取所有特征,共130维,再利用mRMR算法筛选,筛选特征数对比如图5所示。
图5 筛选特征数对比
由图5可知,选择前15个特征分类效果最好,所选特征如表7所示。
表7 所选特征
将所选特征放入RBF_SVM分类器,得到本文所述的完整辅助诊断算法模型,综合各项指标对结果进行评估,评估结果如表8所示。
表8 辅助诊断算法的各项指标 %
如图6和图7所示,10次结果的ACC值约为0.81,说明分类正确的样本数较多,而ROC曲线对应的AUC值在0.89左右,说明分类器性能较优,所以,基于该辅助诊断算法可以实现对草酸钙和无水尿酸结石的准确分类。
图6 辅助诊断算法的ACC,AUC指标
图7 辅助诊断算法的ROC曲线
本文结合临床医学,利用CT图像进行定量分析,并对比机器学习特征提取、选择及分类算法的实验结果,提取结石在二维和三维上的灰度、形状、纹理特征。最终使用mRMR特征选择算法及RBF_SVM分类器,实现患者体内草酸钙和无水尿酸结石的准确分类,ACC和AUC分别达到81.76%和89.03%,比现有结石分析方法更为快速简便,真正实现体内无创,且符合临床统计学结果,为临床医生诊断提供更为有效的参考依据。下一步将优化图像处理和机器学习算法,得到更高的准确率。