陆莉霞,邹俊忠,郭玉成,张 见,王 蓓
1.华东理工大学 控制科学与工程学院,上海200237
2.清影医疗科技(深圳)有限公司,广东 深圳518083
膝关节主要包括韧带、关节面、肌肉、肌腱等多个组织结构,是人体结构中最复杂、最易发生损伤的关节之一。作为临床上一种常见的骨损伤类型,膝关节损伤多发生在半月板、膝关节韧带等位置,患者多表现为不同程度肿胀、疼痛等,且随着现代交通业、建筑业的不断发展,该疾病发生率呈现出逐年增高的趋势,不仅给患者生活带来了极大的不便,还影响到患者的生活质量[1]。
随着X 线、CT、MRI 等影像学检查的普及,膝关节损伤诊断的准确率显著上升,为临床早期治疗提供客观依据。MRI检查具有较高的软组织分辨率,能够采集到更为丰富的图像信息。通过矢状位、横轴位等方位的MRI 影像能够对膝关节软组织及周围结构予以全面的探查,获得详细的、准确的图像,因而MRI检查对前十字韧带撕裂、半月板撕裂的临床诊断表现出较高的准确性[2-3]。近年来,膝关节磁共振成像(MRI)已成为诊断膝关节损伤的首选方法,一个用于分析膝关节MRI影像的自动化系统,可以筛选出高危患者,并帮助临床医生做出更准确的诊断。
传统医学影像分析主要采用边缘检测、纹理特征、形态学滤波以及构建形状模型和模板匹配等方法。Haralick等人[4]使用灰度共生矩阵来提取乳腺组织的纹理特征;Swain 等人[5]使用颜色直方图来提取眼底图像的颜色特征。
MRI影像的多维度、多平面的特性限制了传统医学影像分析方法对膝关节MRI的应用[6],而深度学习方法能够自动学习多层特征,非常适合于辅助诊断医学影像[7]。近年来,深度学习方法已经超越了传统的医学影像分析方法,并在膝关节MRI 影像领域取得了重大进展。Prasoon等人[8]提出了一种新的体素分割系统,该系统集成了三个二维CNN 网络,分别用于MRI 三维影像的xy、yz和zx三个平面的膝关节胫骨软骨的分割。Liu 等人[9]开发了一套基于深度学习的全自动膝关节MRI软骨损伤检测系统,该系统有两个CNN网络组成,第一个CNN 网络进行软骨和骨骼的快速分割;第二个CNN网络评估了分叶软骨组织的结构异常。吴恩达团队[10]首次提出了用于预测膝关节损伤的深度学习模型——MRNet,该模型通过MRI 影像对膝关节常见的损伤进行预测,提高MRI诊断的质量。
本文通过将膝关节矢状面、冠状面和轴向面的MRI影像的多种特征融合,对前交叉韧带撕裂、半月板撕裂和一般异常三种常见的膝关节损伤进行分类预测。通过各个对比实验,验证本文方法的有效性和可行性。
对于膝关节MRI影像的特征提取方法,主要分为传统特征提取方法和深度学习特征提取方法。一般来说,传统方法提取的特征底层语义特征,例如边缘特征和纹理特征。而深度学习如卷积神经网络(Convolution Neural Network,CNN)等,提取的图像特征包含更多的细节信息,能表达高层语义信息。底层特征语义信息与高层语义信息对于膝关节损伤的诊断均有重要意义,因此本文提出用于膝关节诊断的一种多模态特征融合网络,如图1所示。
图1 多模态特征融合网络流程图
首先,提取边缘特征HOG 特征和纹理特征LBP 特征,经contact 融合后利用PCA 选取特征贡献度超过95%的特征作为传统特征;其次,在VGG16网络模型的基础上加入金字塔融合的思想,将多个feature map 的信息融合作为深度特征;最后,通过一个多层神经网络将传统特征和深度特征经进行深度融合,并得到预测概率。
膝关节MRI 影像的一个样本的数据大小为n×256×256,其中n为一组MRI影像序列所包含的静态图片数目。
由于膝关节主要位于MRI的中央位置,因此先将所有静态图片进行裁剪操作,并截取图像中心224×224区域。接着进行归一化操作,像素值被归一化到[0,255]之间,归一化是为了消除其他变换函数对图像变换的影响。在训练阶段,数据进行数据增强,包括水平翻转和旋转角度为25°,水平和垂直方向移动0.1的仿射变换。
1.2.1 HOG特征
HOG特征提取算法是Dalal等[11]提出的一种特征提取算法,核心思想是通过计算图像局部区域梯度,并将每个局部区域中各像素点梯度的方向直方图级联。图2为HOG特征提取算法的基本流程图。
图2 HOG特征提取算法的基本程图
HOG特征提取的具体步骤如下:
(1)为消除光照和噪声的影像,首先对图像进行图像灰度化和伽马标准化处理。
(2)用一维中心对称算子k=[-1 01]及其转置计算图像横坐标和纵坐标方向的梯度,表达式为:
其中,H(x,y)表示像素在(x,y)处的灰度值。得到像素点的梯度幅值和梯度方向为:
(3)将图像分割为若干个小块(block),每个小块由4 个单元(cell)组成,每个单元由8×8 像素组成,块的滑动步长为1个单元。将θ(x,y)在[0,π]区间上分为9个区间,即每20°为一个方向(bin)。单元(cell)中的每一个像素点都为直方图通道进行加权投票,权重为g(x,y),这样就得到每个单元内的9个方向的梯度直方图。
最后按顺序级联所有block 的直方图向量,得到图像的HOG特征μHOG,维度为7 056维。
1.2.2 LBP特征
原始的LBP模式算子为3×3的局部邻域,以邻域中心点像素的灰度值I(c)为基准,对周围其他8 个像素点做二值化处理,比I(c)大的邻域像素值置为1,否则为0。将这8个二进制数按顺序排列得到一个二进制序列即为中心点像素的LBP值。
为对原始的LBP模式进行降维,并在数据量减少的情况下能最好地表示图像的信息,本文采用均匀局部二值模式(Uniform Local Binary Pattern,ULBP)[12]。当某个LBP 所对应的循环二进制数从0 到1 或从1 到0 最多有两次跳变时,该LBP所对应的二进制就称为一个等价模式,其余情况为非均匀模式。通过这样的改进,模式数量由原来的2p种减少为p(p-1)+2,二进制模式的种类大大减少,而不会丢失任何信息。
本文用局部二进制编码直方图表示LBP 图像的特征,将LBP 图像分为8×8 块区域,获取每个区域的LBP编码直方图,继而得到整幅图像的LBP 编码直方图,一张LBP图像共有3 776维特征,记为μLBP。
1.2.3 HOG特征与LBP特征融合
经上述方法提取到的HOG特征和LBP特征均存在特征维度过高,维度之间存在耦合的问题,本文通过主成分分析(PCA)[13]来达到降维的效果。PCA 不仅减少了特征集的维数,同时避免了重要信息的丢失,保持特征集中对方差贡献最大的特征。
设向量X=(x1,x2,…,xn)T,PCA的具体步骤如下:
将HOG 特征和LBP 特征归一化后,经contact 融合后的特征为[μHOG,μLBP]。经PCA 降维后,将特征向量按对应特征值大小从上到下按行排列,前1 325 项的累计贡献率超过95%,为95.43%。因此将原特征集降维到1 325维得到传统特征集μ=(μ1,μ2,…,μ1325)。
近年来,卷积神经网络在医疗图像诊断领域的应用得到了广泛而迅速地发展,并已取得不错的成果。卷积神经网络通过对图像的卷积操作,可自动完成特征的提取,且这些特征拥有更高级的语义信息,鲁棒性更强[14]。
VGG16网络[15]是2014年由牛津大学的视觉几何组(Visual Geometry Group)提出,该网络共有五个卷积块,每个卷积块后都是最大池化层,最后三层为全连接层,VGG16的结构配置见表1。
表1 VGG16结构配置
VGG16 网络利用最后一层特征来进行分类,这种方法的优点是速度快、需要内存少。缺点是仅仅关注深层网络中最后一层的特征,却忽略了其他层的特征,但多种信息的融合可以在一定程度上提升诊断的准确率。由于卷积神经网络中浅层可以学习到比较通用的低级特征,如基本的灰度和边缘信息等,深层可以学到特征的更多细节。因此本文在VGG16 网络的基础上加入Top-down 结构,将提取的特征层逐步搭建为特征金字塔,并进行像素叠加,网络结构如图3所示。
通常深度学习被认为在大量有标注的图像上进行训练才可以达到良好的效果,在医学图像中很难有大量带有标注的病理图像用于分割训练,因此本文采用迁移学习的方法,直接使用预训练好的VGG16 的第二至第五个卷积块卷积块(Block)。将每个卷积块的最后一个卷积层得到的特征映射提取出来,分别为F(n),其中n=2,3,4,5。F(n)经上采样(upsampling)与F(n-1)经1×1 卷积处理后的结果进行像素点融合,再采用3×3 的卷积核对融合结果进行修正,消除了上采样的混叠效应,得到新的特征映射F′(n-1)。金字塔融合的公式为:
获得最后一层融合层的F′(2),依次通过BN 层、ReLU层、自适应最大池化层和全连接层。BN 层不仅能够加快模型的收敛速度,还能提升分类效果[16]经过BN(Batch Normalization)层。设x(k)为F′(2)第k维特征,BN层为x(k)引入了两个可学习的参数γ(k)和β(k),并利用无差估计输出第k维特征:
ReLU 层加入非线性因素的,增加模型的表达能力不够。ReLU激活函数的表达式为:
自适应最大池化层与标准的Max Pooling 区别在于,自适应最大池化层Adaptive Max Pooling 会根据input_size(In)来控制输出output_size(Out)。stride和κernelsize分别为:
全连接层可以看做h×w的全尺寸卷积,h和w为前一层输出的大小,最终得到1 024 维通过卷积神经网络提取的深度特征g=(g1,g2,…,g1024)。
由于提取到的传统特征和深度特征维度均较高,若直接级联两个特征向量作为最后的特征会大大增加时间成本和计算成本,因而在多模态特征融合时要减少融合后特征的维度。同时考虑到不同模态的特征在融合结果中所起到的作用,以及保留多模态的相关性,本节构建了一个多模态特征融合模块,该模块包含一个神经元个数小于特征维数的隐藏层和一个Sigmoid 层,通过最大化特征层的能量比重来训练整个网络,如图4 所示。隐藏层对融合后特征进行非线性降维,Sigmoid 层将融合后特征映射到(0,1)区间,即预测概率。令特征向量x=(μ,g),前向传播公式为:
图3 加入金字塔融合的VGG16网络结构图
图4 多模态特征融合网络图
其中,T为样本数量。当函数ψ(θ)取得最大值时,特征层的能量比重大,隐层的能量小。当在网络内传输数据时,数据流的方向也是能量衰减的方向,多次迭代之后,网络能量呈衰减趋势,网络趋于有序或者概率分布趋于集中。
本文将偏差βi、αj初始值设为服从平均值为0,标准差为0.1的正态分布的随机数,权重γij初始值设为服从平均值为0,标准差为0.01的正态分布很小的随机数,隐藏层神经元设为1 000个。
实验时,所使用的CPU 型号为3.4 GHz、16 核AMD1950x,内存容量为64 GB,GPU 为12 GHz 显存NVIDIA TITAN-V。软件方面,使用Python作为设计语言,配合opencv等视觉库进行代码编写。深度学习模型框架使用Pytorch以完成对模型的构建,训练与预测。
本文网络模型使用Adam优化器,增加模型收敛稳定性和收敛速度,同时也能得到优良的分类结果。学习率经过调试设置为1E-5,权重衰减为0.01,模型运行50轮。
实验过程中使用了公开的膝关节MRI 竞赛数据集[10]对本研究中的网络模型进行训练和验证。本研究的核心任务为检阳,因而本文采用了以下评价指标:准确率Accuracy、Recall 和AUC(Area Under Curve)值定量的评估不同模型的性能。
准确率为模型预测正确的样本占所有参与预测的样本的比例。定义为:
其中,TP为真阳性样本数量,TN为真阴性样本数量,FP为假阳性样本数量,FN为假阴性样本数量。
AUC 值被定义为ROC 曲线下与坐标轴围成的面积。由于ROC曲线一般都处于y=x这条直线的上方,所以AUC 的取值范围在0.5 和1 之间。AUC 越接近1.0,检测方法真实性越高。
本次研究中的膝关节MRI 影像数据集来自于斯坦福大学医学中心整理的从2001 年1 月1 日至2012 年12 月31 日期间内的1 370 个膝关节核磁共振检查数据集——MRNet数据集[18]。该数据集包含1 104例(80.6%)异常(Abnormal)检查,其中前交叉韧带撕裂(ACL tear)319例(23.3%),半月板撕裂(Meniscus tear)508例(37.1%)。194例(38.2%)患者同时发生前交叉韧带撕裂和半月板撕裂。每个样本包括矢状面加权序列(sagittal plane)、冠状面加权序列(coronal plane)和轴向面加权序列(axial plane)。这些序列的图像数量从17 幅到61 幅不等(平均31.48幅,SD7.97幅)。
数据集被分为训练集(1 130个样本,1 088名患者)和验证集(120个样本,111名患者),训练集与测试集的数据分布见表2。
表2 训练集与测试集的数据分布
若将矢状面加权序列、冠状面加权序列和轴向面加权序列同时放入上述多模态特征融合网络中,不仅对硬件要求较高,运行时间较长,而且不同序列的特征表达差异较大,预测准确度不高。
因而本研究针对每个任务(异常、前交叉韧带撕裂、半月板撕裂)和序列类型(矢状面、冠状面、轴向面)训练了不同的9个模型。
在验证阶段,采用回归模型得到异常、前交叉韧带撕裂、半月板撕裂的预测概率,回归模型结构图如图5所示。
图5 回归模型结构图
将同一个样本的矢状面、冠状面和轴向面加权序列同时放入上述多模态特征融合网络中,分别得到在该序列下,样本为异常样本的概率P(A)、P(B)和P(C)。然后,将P(A)、P(B)和P(C)分别乘以一个回归系数,相加后通过Sigmoid 函数得到在三种序列下,此样本为异常样本的预测概率,计算公式为:
为了寻找最佳参数wA、wB和wC,本文采用梯度上升的方法,将梯度算子移动的步长记为α,则梯度算法的迭代公式为:
同样地,可以获得此样本为前交叉韧带撕裂样本和半月板撕裂样本的预测概率。
特征提取不仅仅是分类和预测的前提,也是分类准确率的保证,而高效准确的特征可以提高后续分类和预测的精确度。本文数据经相同的预处理后,每种膝关节MRI序列提取的传统特征(Model a)、CNN特征(Model b)和多模态融合特征在每个任务下的实验结果见表3。
由表3 中Model a 的结果可知,仅仅依靠MRI 影像的纹理和边缘特征,在膝关节损伤预测任务中存在一定的局限性,尤其在半月板撕裂的预测上,最高的准确率仅有69.17%,无法很好地辅助临床医生进行预测。而与传统方法相比,深度学习方法在膝关节MRI影像诊断领域具有明显优势。由表3中Model a和Model b的结果可知,深度学习方法的准确率、召回率和AUC值均有大幅提升,其中在异常预测任务中召回率最高已达到97.89%,同时在传统方法有所欠缺的半月板撕裂预测任务中,深度学习方法的准确率均高于76%,充分显示了深度学习方法能够更好地诊断出阳性的样本,同时AUC值也更高,因而在偏态的样本中更稳健。
表3 各模型效果对比
本文模型在卷积神经网络中融入低级语义信息,丰富了模型的特征表达。由表3中Ours的结果可知,在融合了传统特征与深度学习特征后,各性能指标均得到一定程度的提升,尤其是准确率和AUC值。在异常、前交叉韧带撕裂、半月板撕裂这三项任务中,最高准确率为92.50%、92.50%和81.83%,最高AUC 值为0.961 3、0.974 2 和0.857 4。同时,本文模型对一般异常的阳性样本的检出率均高于94%,对前交叉韧带撕裂和半月板撕裂阳性样本的检出率也较高。这对于优先发现高危患者并协助临床医生进行诊断有着重要意义。
由表3 中Ours 的结果可知,在一般异常预测任务中,轴向面序列表现最佳,准确率、召回率和AUC 值分别为92.50%、97.95%和0.961 3;在前交叉韧带撕裂预测任务中,矢状面序列更有优势,准确率、召回率和AUC值分别为92.50%、94.44%和0.972 8;在半月板撕裂预测任务中,轴向面序列表现较好,准确率、召回率和AUC值分别为81.83%、82.92%和0.846 5。此结果与半月板位于两侧胫骨与股骨外侧髁间[19]而前交叉韧带位于膝关节滑膜外侧[20]的膝关节结构相符合,因此,本文实验结果具有一定的科学性。
鉴于每个任务均只有阳性和阴性两种可能,为了更好地展示回归模型的效果,将回归模型预测概率大于或等于0.5 的认为是阳性样本,小于0.5 的认为是阴性样本,则在一般异常、前交叉韧带撕裂和半月板撕裂下的实验结果如表4。
表4 回归模型各性能指标
由表4 可知,本文模型在一般异常、前交叉韧带撕裂和半月板撕裂预测任务中均有较好的性能指标。尤其是一般异常和前交叉韧带撕裂预测任务,准确率和召回率均超过90%,AUC 值均高于0.94,说明本次研究所应用的方法能够在这两项膝关节损伤预测中为临床医生提供可靠的基础诊断。由于半月板主要由纤维软骨构成,而MRI影像对软骨组成成分具有局限性[21]。在这种情况下,本文模型在半月板撕裂预测任务中的AUC值仍达到了0.847 9。综合来看,本文方法在膝关节损伤诊断方面具有一定的价值。
表5 不同膝关节损伤检测模型性能比较
表5 对比了本文算法与MRNet 竞赛[18]前三名结果和现有文献中最好结果的平均AUC值。
可以发现,在对膝关节损伤进行预测时,表5中6种深度学习网络AUC 值均高于90%,都能够取得较好的预测结果,验证了深度学习强大的特征学习能力。Bien等[10]在论文中提出的MRnet 以及mrnet-baseline 都是以alexnet 模型为基础的,AUC 值虽已达到0.917,但受alexnet模型本身深度和广度的限制[22],卷积层能学习到的特征的数量和种类有限。本文模型以更深层的vgg16模型为基础,并在该卷积神经网络的基础上加入了多模态的融合特征,包括多层卷积层的融合特征以及MRI影像的纹理和边缘特征,让模型有了更强大的表达能力,以此来获得更好的预测效果。因此本文模型的平均AUC值略高与其他深度学习网络,为0.920。
本文提出了用于膝关节诊断的一种多模态特征融合网络。首先,通过对膝关节MRI影像的预处理,从传统与深度学习两方面提取膝关节损伤的多模态特征;然后,通过多层神经网络对特征进行快速的相关性融合;最后,用logistic 回归整合各个模型的概率。通过不同的实验,对上述所提到的三种模型进行了比较,同时对回归模型进行了分析,结果显示本文方法在膝关节损伤预测方面具有较大的优势和较高的实用价值。在未来的研究工作中,将尝试更多模态的特征融合,例如MRI影像在时间维度上的特征。同时将此应用在膝关节MRI 的算法应用在其他医学图像类型和不同人体器官的诊断上。