基于非下采样轮廓波变换和卷积神经网络的骨龄评估方法与流程

本发明涉及医学图像分析领域，具体涉及一种基于非下采样轮廓波变换和卷积神经网络的骨龄评估方法。

背景技术：

作为儿科放射学中的一种常用技术，骨龄评估主要用于判断儿童骨骼发育的成熟程度，对内分泌紊乱诊断、生长激素治疗监控、青少年最终身高预测等多种临床问题具有十分重要的意义。目前，临床上最为流行的骨龄评估方法是基于受试者左手及手腕的x光图像，医生通过分析x光图像中相应骨骼结构的发育程度，完成对受试者的骨龄评估。然而，骨龄评估结果的准确性严重依赖于医生的经验和水平，不同的医生评估的结果可能差异很大，即使对于同一位医生，其多次评价的结果也可能存在一定的差别。

近年来，基于计算机辅助的自动骨龄评估方法成为医学图像分析领域一个活跃的研究方向，其基本思想是采用机器学习方法，使计算机具有“读”骨龄片的能力，能够自动给出一个骨龄估计值。借助于计算机强大的计算能力，计算机辅助自动骨龄评估可以在很短的时间完成评估任务，很大程度上减轻医生的工作负担，为医生提供参考和帮助。十几年来，多种自动骨龄评估方法被国内外研究学者提出，这些方法通常将骨龄评估问题建模为一个分类问题或者回归问题，通常包括手部区域分割，感兴趣区域提取、特征提取、分类器/回归器学习等步骤。最近，深度学习理论和方法在越来越多的医学图像分析问题领域获得了巨大的成功。在骨龄评估方向，近两年也逐渐开始出现基于深度学习(主要是卷积神经网络)的相关方法，并且展现出相对于传统方法的优势，为解决该问题提供了新途径。可参考文献：h.lee,s.tajmir,j.lee,etal.,“fullyautomateddeeplearningsystemforboneageassessment”,journalofdigitalimaging,vol.30,no.4,pp.427-441,2017和c.spampinato,s.palazzo,d.giordano,etal.,“deeplearningforautomatedskeletalboneageassessmentinx-rayimages”,medicalimageanalysis,vol.36,pp.41-51,2017。

在目前已有的基于深度学习的骨龄评估方法中，神经网络的输入都是空间域x光图像，而没有进行特征的预提取。然而，在众多医学图像相关应用中，带标签(专家标注)的影像获取成本很高，数据量往往非常有限，这给网络的训练带来了很大困难。即使通过迁移学习的方式利用其它自然图像任务中预训练的网络模型，在样本规模很小的情况下也难以保证网络训练的水平，限制了评估系统的精度。

技术实现要素：

本发明为克服现有技术存在的问题，提供了一种基于非下采样轮廓波变换和卷积神经网络的骨龄评估方法，以期能解决现有基于深度学习的方法在小规模数据集上进行网络训练困难的问题，提升网络的泛化性能，从而为临床应用提供更加准确可靠的评估方法。

本发明为解决上述问题采用如下技术方案：

本发明一种基于非下采样轮廓波变换和卷积神经网络的骨龄评估方法的特点在于，包括如下步骤：

步骤1：以手部x光图像作为原始输入图像并进行插值和边界连续型延拓的空间尺寸处理，得到归一化大小的手部x光图像；

步骤2：对所述归一化后的手部x光图像进行非下采样轮廓波变换处理，得到低频分量和不同尺度下的高频方向子带；

步骤2.1：对所述归一化后的手部x光图像进行非下采样金字塔分解，得到低频分量c0和l个尺度的高频分量{c1,c2,…,ci,…,cl}，其中ci表示第i个尺度的高频分量，1≤i≤l，l为分解层数；

步骤2.2：对第i个尺度的高频分量ci进行非下采样方向滤波，得到第i个尺度下的高频方向子带其中ci,j表示第i个尺度下的第j个高频方向子带，1≤j≤ki，ki为第i个尺度下高频方向子带的数目，从而得到l个尺度下的高频方向子带；

步骤3：第i个尺度下的高频方向子带进行层叠处理，得到第i个立体数据di，从而得到l个立体数据{d1,d2,…,di,…,dl}；

将低频分量c0复制并扩展成一个三通道的立体数据，并记为第l+1个立体数据dl+1；

由l个立体数据{d1,d2,…,di,…,dl}和第l+1个立体数据dl+1构成立体数据集合{d1,d2,…,dk,…,dl+1}，1≤k≤l+1；

步骤4：构造一个包含l+1个通道的卷积神经网络，且每个通道包含若干个卷积层和最大池化层；

将第k个立体数据dk对应输入到第k个通道中，得到第k组输出特征图fk＝conv(dk)，其中conv表示卷积网络的前向过网络算子，从而得到l+1组特征图{f1,f2,…,fk,…,fl+1}；

步骤5：将所述l+1组特征图{f1,f2,…,fk,…,fl+1}进行层叠处理，得到特征图集合f；

步骤6：构造包含若干个全连接层的回归网络，且令最后一层的神经元数量为1；

定义回归网络的前向过网络算子为reg；

将所述特征图集合f输入到所述回归网络中，得到骨龄预测结果r＝reg(f)。

与现有基于深度学习的骨龄评估方法相比，本发明方法对输入x光骨龄图像进行非下采样轮廓波变换，获得图像在不同尺度下的特征表示，将它们各自输入到一个卷积网络进行进一步的特征提取，最后将各个尺度下卷积网络输出的特征融合在一起用于骨龄预测，其特色之处在于：

第一，使用非下采样轮廓波变换可以得到原始图像一组丰富的特征，是一种特征预提取的方式，同时也可视为对原始数据集的一种扩充方法，有效克服了现有深度学习方法在小规模数据集下进行网络训练面临的困难，提升了网络的泛化性能。

第二，使用非下采样轮廓波变换可以对原始空间域图像的特征信息实现很好的分离，不同尺度的特征可以分开处理，再在后端进行融合，这种方式显著增加了卷积神经网络结构设计的灵活性，提高了特征提取的质量，进而提升了骨龄评估系统的性能。

附图说明

图1为本发明骨龄评估方法的流程图；

图2为本发明一幅x光骨龄图像示意图；

图3为本发明x光图像尺寸归一化结果示意图；

图4为现有技术中非下采样轮廓波变换原理示意图；

图5为本发明x光图像进行非下采样轮廓波变换的结果示意图；

图6为本发明使用的vggnet-16网络卷积部分结构示意图；

图7为本发明骨龄评估网络结构示意图。

具体实施方式

图1为本发明一种基于非下采样轮廓波变换和卷积神经网络的骨龄评估方法流程图，共包括以下步骤：

步骤1：以手部x光图像作为原始输入图像并进行插值和边界连续型延拓的空间尺寸处理，得到归一化大小的手部x光图像；

在本实施例中，以美国南加州大学提供的公开骨龄评估数据集digitalhandatlas中的一幅骨龄图像为例，如图2所示，受试者是一位实际年龄为10.49岁的亚裔女孩，该图像的原始空间分辨率是1526*1893像素。由于后续卷积网络需要输入图像的空间分辨率固定，同时考虑到计算效率和存储代价，故需要对原始输入图像的尺寸进行归一化。这里，将原始输入图像归一化到256×256像素，同时通过连续型边界延拓保持图像内容形状不变。具体做法为，将图像先进行保持长宽比的插值，使长边为256像素，然后对短边两端进行连续型边界延拓至256像素。图3为图2中原始x光图像进行尺寸归一化后的结果。

步骤2：对归一化后的手部x光图像进行非下采样轮廓波变换处理，得到低频分量和不同尺度下的高频方向子带；

步骤2.1：对归一化后的手部x光图像进行非下采样金字塔分解，得到低频分量c0和l个尺度的高频分量{c1,c2,…,ci,…,cl}，其中ci表示第i个尺度的高频分量，1≤i≤l，l为分解层数；

步骤2.2：对第i个尺度的高频分量ci进行非下采样方向滤波，得到第i个尺度下的高频方向子带其中ci,j表示第i个尺度下的第j个高频方向子带，1≤j≤ki，ki为第i个尺度下高频方向子带的数目，从而得到l个尺度下的高频方向子带；

非下采样轮廓波变换(non-subsampledcontourlettransfrom,nsct)是一种代表性的图像多尺度几何变换方法，具有多尺度、多方向、平移不变等优良特性。非下采样轮廓波变换的实现主要基于非下采样金字塔(non-subsampledpyramid,nsp)分解和非下采样方向滤波器组(non-subsampleddirectionalfilterbank，nsdfb)分解两部分。nsp分解用于获得输入图像的多尺度表示，第一层nsp分解作用于原始输入图像，得到和输入图像尺寸相同的一个低频系数图和一个高频系数图，而后通过对上一层的低频系数图迭代地进行nsp分解，得到图像由细到粗的的多尺度分解结果。因此，如果进行l层nsp分解，可以得到l幅高频系数图和1幅低频系数图。对每一个尺度的高频系数图，再使用nsdfb分解可以获得相应尺度的多方向表示，每个方向得到的系数图和原始图像的大小也相同。通常情况下，每一层的方向数是2的幂次，且越细的尺度分解的方向数越多。图4为非下采样轮廓波变换原理示意图。

在本实施例中，对输入图像进行四层nsp分解，即l＝4，对四个尺度的高频分量(由细到粗)分别进行16、8、8、4方向的nsdfb分解。图5为图3中的x光图像进行上述非下采样轮廓波变换的结果示意图。由图5可以看出，不同尺度的系数图包含的特征信息存在着显著的区别。

步骤3：第i个尺度下的高频方向子带进行层叠处理，得到第i个立体数据di，从而得到l个立体数据{d1,d2,…,di,…,dl}；

将低频分量c0复制并扩展成一个三通道的立体数据，并记为第l+1个立体数据dl+1；

由l个立体数据{d1,d2,…,di,…,dl}和第l+1个立体数据dl+1构成立体数据集合{d1,d2,…,dk,…,dl+1}，1≤k≤l+1；

在本实施例中，四个尺度高频数据的高度维度分别为16、8、8、4。

步骤4：构造一个包含l+1个通道的卷积神经网络，且每个通道包含若干个卷积层和最大池化层；

将第k个立体数据dk对应输入到第k个通道中，得到第k组输出特征图fk＝conv(dk)，其中conv表示卷积网络的前向过网络算子，从而得到l+1组特征图{f1,f2,…,fk,…,fl+1}；

在本实施例中，采用视觉识别领域非常流行的vggnet-16网络的卷积部分作为每个尺度系数图的特征提取模型。图6为vggnet-16网络卷积部分结构示意图，网络输入大小为224×224像素，共包含13个卷积层和5个最大池化层，所有卷积层的卷积核空间大小均为3×3，统一用“conv3”表示，其后的数字表示相应层输出的特征图数量，由64逐渐增加到512，最大池化层用“max-pool”表示。

步骤5：将l+1组特征图{f1,f2,…,fk,…,fl+1}进行层叠处理，得到特征图集合f；

在本实施例中，卷积网络共包含5个通道，每个通道输出512幅特征图，故合并后共包含2560幅特征图。

步骤6：构造包含若干个全连接层的回归网络，且令最后一层的神经元数量为1；

定义回归网络的前向过网络算子为reg；

将特征图集合f输入到回归网络中，得到骨龄预测结果r＝reg(f)。

在本实施例中，使用一个包含三个全连接层的网络结构作为回归网络，这三个全连接层依次包括1024、128和1个神经元，最后一个神经元的输出结果即为骨龄预测值。

上述步骤以一个端到端(end-to-end)的卷积神经网络结构实现，网络结构示意图如图7所示，输入为待评估x光图像，输出为骨龄预测结果。在训练过程中，使用专家标注的数据作为训练样本集，并进行多视窗、水平镜像等数据扩充操作，将真实值与预测值之间的误差平方作为损失函数。具体地，令θ表示网络所有连接权重的参数集合，损失函数定义为

其中，n表示一批(batch)训练样本包含的样本数量，xi表示第i个训练样本，f(xi,θ)表示网络预测的骨龄值，yi表示第i个训练样本的真实骨龄值，由专家提供。利用误差反向传播机制实现网络训练，使用随机梯度下降法迭代计算得到所有连接的权重。在测试过程中，使用中心视窗或多视窗融合的方式，将测试样本输入到训练好的网络中，得到骨龄预测值。

在本实施例中，使用美国南加州大学提供的公开骨龄评估数据集digitalhandatlas进行网络训练和性能验证。该数据集共包含近1400幅0-18岁受试者的手部x光骨龄图像，覆盖不同的性别和种族。对于每幅图像，两名儿科专家各自给出一个评估结果。本实施例的相关设置如下：将图像数据均分为五份，采用留一法进行实验，其中四份用于训练，一份用于测试，循环五次进行交叉验证；在caffemodelzoo上下载在imagenet数据集上预训练的vggnet-16网络模型作为预训练模型，在此基础上进行微调(由于网络输入数据不再是三通道的rgb图像，第一个卷积层的参数使用随机初始化的方式)；将256×256大小的输入数据进行9个视窗、大小为224×224(vggnet-16网络输入大小)的采样和水平镜像以扩充数据集，防止过拟合；随机梯度下降法的动量因子设为0.9，权重衰减系数设为0.0005；网络训练共迭代80个周期(epoch)，前60个epoch的学习速率为0.0001，后20个epoch的学习速率设为0.00001；对每个样本，两名专家评估结果的平均值作为真实值用于网络训练；在测试时，将数据的中心窗口输入到训练好的网络中，得到骨龄预测值；在性能验证时，将网络预测结果与每名专家的评估结果分别比较，计算平均绝对值误差作为性能指标。

为了体现本发明所提供方法的有效性，这里将上述实施例提供的方法与三种现有的基于深度学习的骨龄评估方法在digitalhandatlas数据集上进行性能比较。表1列出了本发明实施例提供的方法与三种现有深度学习骨龄评估方法的性能测试结果。三种对比方法均以原始空间域图像作为网络输入，而没有进行特征预提取，区别在于使用的网络结构不同，分别为vggnet-16、googlenet和bonet，相关实验设置和本发明实施例提供的方法一致。由表1可以看出，对于任何一名评审专家的评估结果，本发明实施例提供的方法都能给出比其他三种方法更为接近的预测值，体现出本发明方法的有效性。此外，本发明方法具有较高的计算效率，在一台配有4.0ghzcpu和nvidiagtx1080gpu的计算机上，评估一幅x光骨龄图像的时间花费约为0.03秒，可以达到实时性要求。

表1本发明实施例提供的方法与三种现有深度学习骨龄评估方法的性能比较结果

(误差单位：年)