技术领域
[0001] 本
发明涉及医学影像处理技术领域，具体涉及一种基于可形变配准和DCNN的宫颈MRI图像分割方法。
背景技术
[0002] 在失去手术机会的情况下，
放射治疗对于
宫颈癌是至关重要的。与CT(Computed Tomography,计算机
断层扫描)相比，MRI(Magnetic Resonance Imaging,
磁共振成像)不仅没有
电离辐射还可以提供较好的软组织
对比度，越来越多地被用于
放射治疗。精确
定位宫颈及其周围结构是MRI引导宫颈癌放射治疗的必要条件。手工分割多个组织器官十分耗时，而且可能会存在不同观察者的差异。由于具有灰度不均匀性且器官的解剖变异较大，宫颈MRI图像的计算机辅助分割任务仍是一项具有挑战性的工作。如图1所示，该挑战主要与灰度不均匀性和器官形状变化有关。
[0003] 为了解决这些问题，通常采用atlas-to-target(模板图像到目标图像)的可形变配准方法。Staring等在文献中采用KNNG(k-Nearest Neighbors Graph,k-最近邻图)的方法来实现α-MI(Multi-feature Mutual Information,多特征互信息)，对宫颈MRI图像进行配准。与MI(Mutual Information,标准的互信息)方法相比，他们使用笛卡尔特征集的方法明显提高了配准
精度。Lu等在文献中将自动分割、非刚性配准和
肿瘤检测集成到贝叶斯
框架中，该框架需要使用
水平集表示来构造形状先验。Berendsen等在患者间的
图像配准中，将统计形状模型视为自由形变框架内的一个惩罚项。然而，由于
训练数据的形状和大小变化较大，建立统计形状模型并不是一项简单的任务。
[0004] 因此，亟需一种精度更高的实现宫颈MRI图像分割的技术方案。
发明内容
[0005] 针对
现有技术中存在的
缺陷，本发明的目的在于提供一种基于可形变配准和DCNN的宫颈MRI图像分割方法，以获得更高的分割精度。
[0006] 为达到以上目的，本发明采取的技术方案是：一种基于可形变配准和DCNN的宫颈MRI图像分割方法，包括以下步骤：
[0007] 采用目标图像的自相似性加权初始生成图的边，计算基于最小距离树的加权α-MI度量；
[0008] 采用加权α-MI度量进行目标图像和模板图像之间的可形变配准，期间采用
梯度下降法实现参数寻优，通过配准结果对模板标签图像进行形变，获得粗分割结果；
[0009] 采用粗分割结果的符号距离场定位潜在的边界区域，采用经过迁移学习的DCNN模型进行所述潜在边界区域的分类预测；
[0010] 对DCNN模型输出的预测结果进行边界细化提炼，得到精分割结果。
[0011] 在上述技术方案的
基础上，计算基于最小距离树的α-MI度量，具体包括以下步骤：
[0012] 假设z(xi)＝[z1(xi)，…，zd(xi)]表示一个包含了点xi所有特征的d维度向量，Tμ(·)是参数为μ的形变模型，zf(xi)为目标图像上点xi的
特征向量，zm(Tμ(xi))为坐标转换后的点Tμ(xi)在模板图像上的特征向量，zfm(xi，Tμ(xi))为两个特征向量的拼接：[zf(xi)，zm(Tμ(xi))]，依据下面的公式构造三个最小距离树：
[0013]
[0014]
[0015]
[0016] 其中，||·||表示欧氏距离，γ∈(0，d)，n表示
采样点的个数，i和j是两个不同点xi和xj的下标；
[0017] 基于最小距离树的α-MI度量的计算公式为：
[0018]
[0019] 其中α＝(d-γ)/d，Lf,Lm和Lfm分别为上述构造的三个最小距离树的长度。
[0020] 在上述技术方案的基础上，采用目标图像的自相似性加权初始生成图的边，计算基于最小距离树的加权α-MI度量，具体包括以下步骤：
[0021] 构造一个直方图描述子，将以当前
体素为中心的图像
块投影到二维直方图中，该直方图的坐标轴是中心体素到样本点的径向距离和样本点的强度值；
[0022] 采用Moran′I空间自相关系数将自相似估计限制在结构信息丰富的图像区域，用于计算该系数的图像块大小与用于构造直方图描述子的图像块大小相同；
[0023] 选择Moran′I系数值大于整体标准差的
像素点，使用两个直方图描述子之间的陆地移动距离计算两个像素的相似性ωij；
[0024] 使用自相似性度量ωij对初始生成图的边进行加权，改良公式(3)为：
[0025]
[0026] 其中zfm(xi，Tμ(xi))为两个特征向量的拼接，n表示采样点的个数，||·||表示欧氏距离，γ∈(0，d)；
[0027] 将公式(5)替代公式(3)代入到公式(4)中即可获得基于最小距离树的加权α-MI度量。
[0028] 在上述技术方案的基础上，采用加权的α-MI度量进行目标图像和模板图像之间的可形变配准，期间采用梯度下降法实现参数寻优，通过配准结果对模板标签图像进行形变，获得粗分割结果，具体包括以下步骤：
[0029] 采用梯度下降法实现参数寻优，需要计算加权的α-MI度量相对于形变模型参数μ的梯度解析表达式：
[0030]
[0031] 其中α＝(d-γ)/d，γ∈(0，d)，d是特征向量的维度，Lm和Lfm为公式(2)和(3)定义的两个最小距离树；
[0032] 为了更简洁的表达，预先定义 代表两点特征向量的距离：
[0033]
[0034] 其中zfm(xi，Tμ(xi))为两个特征向量的拼接，从公式(5)可以推导出Lfm的梯度为：
[0035]
[0036] 进一步地有，
[0037]
[0038] 其中， 是模板特征图像的空间导数， 是该变换的雅可比矩阵；
[0039] 采用类似的方法得到Lm的梯度与公式(9)一起代入到公式(6)中进行计算；
[0040] 梯度下降法通过
迭代的方式获得最佳的形变场 计算公式有，
[0041] μk+1＝μk+ak·dk (10)
[0042] 其中ak为步长因子，dk为通过公式(6)计算的加权α-MI度量的梯度，μk为第k步的形变模型参数；
[0043] 采用最佳形变场 对模板标签图像进行形变，获得目标图像的粗分割结果。
[0044] 在上述技术方案的基础上，采用粗分割结果的符号距离场定位潜在的边界区域，采用经过迁移学习的DCNN模型进行该潜在边界区域的分类预测，具体包括以下步骤：
[0045] 在预处理过程中，每幅图像的灰度值被归一化为[0,255]；从轴向视
角对每个二维切片实施精分割；
[0046] 选择VGG-19网络执行从自然图像到MRI图像的迁移学习，采用ImageNet
数据库中的自然图像对VGG-19模型进行预训练；
[0047] 从当前任务的训练切片中提取图像块对已预训练好的VGG-19模型进行微调；在每个训练切片中，通过人工手动分割结果，边界区域被定义为符号距离值在-5到5之间的一组像素；在此边界区域中随机选取
种子像素，以每个种子像素为中心提取一个图像块，并重复三次作为训练样本；
[0048] 对于上一阶段的粗分割结果，通过符号距离场定位潜在的边界区域，该边界区域被定义为符号距离值在-5到5之间；对于该区域的每个像素，以它为中心提取图像块，输入到已训练好的VGG-19网络模型中，将该像素预测分类为前景或背景。
[0049] 在上述技术方案的基础上，对VGG-19模型进行预训练，具体包括以下步骤：
[0050] 在拥有1000类自然图像的ImageNet数据库上进行VGG-19网络的预训练；
[0051] VGG-19由16个卷积层、5个
池化层和3个全连接层组成；保留最后一个全连接层的两个神经元，并移除其他输出神经元。
[0052] 在上述技术方案的基础上，对DCNN模型输出的预测结果进行边界细化提炼，得到精分割结果，具体包括以下步骤：
[0053] 采用中值滤波以去除单独的刺和洞；
[0054] 移除10％的边界点；
[0055] 用立方的b样条拟合剩下的边界点，形成精分割。
[0056] 与现有技术相比，本发明的优点在于：
[0057] 本发明提供了一种可形变配准和
深度神经网络相结合的由粗到精的分割策略。首先对采用最小距离树的多特征互信息实施改良，利用这种新型度量完成可形变配准实现对潜在边界区域的粗分割，这种改良策略是对初始生成图的边进行自相似估计加权，使其更加紧凑，利于配准参数寻优；其次，通过迁移学习的方法将在自然图像上性能良好的像素分类器深度神经网络应用到当前任务中，从而重新预测潜在边界区域中每个像素的类别标签，完成精分割，在宫颈MRI图像上提高了分割性能，获得了更高的分割精度。
附图说明
[0058] 图1为背景技术中宫颈MRI图像示例图；
[0059] 图2为本发明
实施例的基于可形变配准和DCNN的宫颈MRI图像分割方法的流程示意图；
[0060] 图3为本发明实施例的定位边界区域取图像块的示意图；
[0061] 图4为本发明实施例的不同方法在不同解剖结构下的配准精度DSC箱线图，*表示该列与前一列相比有统计显著性差异；
[0062] 图5为本发明实施例的分割结果对比展示图。
具体实施方式
[0063] 以下结合附图对本发明的实施例作进一步详细说明。
[0064] 参见图2所示，本发明实施例提供一种基于可形变配准和DCNN的宫颈MRI图像分割方法，包括以下步骤：
[0065] S1、采用目标图像的自相似性加权初始生成图的边，计算基于最小距离树的加权α-MI度量；
[0066] S2、采用加权α-MI度量进行目标图像和模板图像之间的可形变配准，期间采用梯度下降法实现参数寻优，通过配准结果对模板标签图像进行形变，获得粗分割结果；
[0067] S3、采用粗分割结果的符号距离场定位潜在的边界区域，采用经过迁移学习的DCNN模型进行所述潜在边界区域的分类预测；
[0068] S4、对DCNN模型输出的预测结果进行边界细化提炼，得到精分割结果。
[0069] S1、基于MST(Minimal Spanning Tree,最小距离树)的加权α-MI
[0070] 假设z(xi)＝[z1(xi)，...，zd(xi)]表示一个包含了点xi所有特征的d维度向量，Tμ(·)是参数为μ的形变模型，zf(xi)为目标图像上点xi的特征向量，zm(Tμ(xi))为坐标转换后的点Tμ(xi)在模板图像上的特征向量，zfm(xi，Tμ(xi))为两个特征向量的拼接：[zf(xi)，zm(Tμ(xi))]。
[0071] 三个MST图可以依据下面公式进行构造：
[0072]
[0073]
[0074]
[0075] 其中，||·||表示欧氏距离，γ∈(0，d)，n表示采样点的个数，i和j是两个不同点xi和xj的下标；
[0076] 基于MST的α-MI度量(称为aMI-MST)可以表达为：
[0077]
[0078] 其中α＝(d-γ)/d，Lf，Lm和Lfm分别为上述构造的三个MST图的长度。
[0079] 为了解决局部极值问题，利用目标图像的自相似性来加权初始生成图的边。首先，构造一个直方图描述子，将以当前体素为中心的图像块投影到二维直方图中，该直方图的坐标轴是中心体素到样本点的径向距离和样本点的强度值。然后，利用Moran′I空间自相关系数将自相似估计限制在结构信息丰富的图像区域。用于计算该系数的图像块大小(r＝4)与用于构造直方图描述子的图像块大小相同。对于在一个局部邻域内的自相似估计，选择那些Moran′I系数值大于整体标准差的像素点。最后，使用两个直方图描述子之间的EMD(Earth Mover’s Distance,陆地移动距离)来计算两个像素的相似性ωit。
[0080] 使用自相似性度量ωij对初始生成图的边进行加权，改良公式(3)为：
[0081]
[0082] 其中zfm(xi，Tμ(xi))为两个特征向量的拼接，n表示采样点的个数，||·||表示欧氏距离，γ∈(0，d)；
[0083] 将公式(5)替代公式(3)代入到公式(4)中即可获得基于MST的加权α-MI度量(称为aMI-WMST)。
[0084] S2、使用可形变配准进行粗分割
[0085] 在第一阶段，粗分割使用了一种atlas-to-target的配准方法进行。通常，atlas由模板灰度图像和模板标签图像组成。模板标签图像通常是由专家手动分割模板灰度图像而得。采用加权的α-MI度量进行目标图像和模板图像之间的可形变配准，期间采用梯度下降法实现参数寻优。通过配准结果对模板标签图像进行形变，可以获得目标解剖结构的初始轮廓。同时为α-MI度量选择一组笛卡尔图像结构特征，其中包含了一阶和二阶梯度。无论是旋转还是平移，这些特征都是不变的。这些不变量可以用范围为σ＝1,2的高斯梯度计算。
[0086] 当随机梯度下降策略被应用于图像配准的参数优化问题时，加权的α-MI度量相对于变换参数μ的导数是必需的。因此，有：
[0087]
[0088] 其中α＝(d-γ)/d，γ∈(0，d)，d是特征向量的维度，Lm和Lfm为公式(2)和(3)定义的两个MST图；
[0089] 为了更简洁的表达，预先定义 代表两点特征向量的距离；
[0090]
[0091] 其中zfm(xi，Tμ(xi))为两个特征向量的拼接，从公式(5)可以推导出Lfm的梯度为：
[0092]
[0093]
[0094] 进一步地有，
[0095]
[0096] 其中， 是模板特征图像的空间导数， 是该变换的雅可比矩阵；
[0097] 采用类似的方法得到Lm的梯度与公式(9)一起代入到公式(6)中进行计算；
[0098] 梯度下降法通过迭代的方式获得最佳的形变场 计算公式有，
[0099] μk+1＝μk+ak·dk (10)
[0100] 其中ak为步长因子，dk为通过公式(6)计算的加权α-MI度量的梯度，μk为第k步的形变模型参数。
[0101] 采用最佳形变场 对模板标签图像进行形变，获得目标图像的粗分割结果。
[0102] S3、使用深度
卷积神经网络进行精分割
[0103] 在预处理过程中，每幅图像的灰度值被归一化为[0,255]。由于切片间的厚度比xoy平面像素的
分辨率大得多，从轴向视角对每个二维切片实施精分割。
[0104] 因为VGG-19网络在大量自然图像上具有良好的表现，所以选择它来执行从自然图像迁移学习到MRI图像的分割任务。VGG-19由16个卷积层、5个池化层和3个全连接层组成。每个卷积层的输入与多个3×3
内核进行卷积从而得到一组特征图。VGG-19网络的预训练在拥有1000类自然图像的ImageNet数据库上完成。为了精分割，保留了最后一个全连接层的两个神经元，并移除其他输出神经元。
[0105] 从训练切片中提取图像块对预训练好的VGG-19模型进行微调。如图3所示，(a)为人工手动的标记结果；(b)为该标记结果的符号距离场；(c)为在对应的灰度图像上提取训练块的操作，在每个训练切片中，通过人工手动分割结果，边界区域被定义为符号距离值在-5到5之间的一组像素。在此边界区域中随机选取种子像素，以每个种子像素为中心提取一个48×48的图像块，并重复三次作为训练样本。将种子像素的所属类别视为本次训练样本的类标签，训练的图像块数量约为100,000个，其中正样本和负样本各占一半。
[0106] 对于上一阶段的粗分割结果，通过符号距离场定位潜在的边界区域，该边界区域同样被定义为符号距离值在-5到5之间；对于该区域的每个像素，以它为中心提取图像块，输入到已训练好的VGG-19网络模型中，将该像素预测分类为前景或背景。
[0107] S4、对DCNN模型输出的预测结果进行边界细化提炼
[0108] 由于VGG-19模型输出的分割结果可能包含不规则区域，因此需要进一步的边界提炼。在此过程中，首先采用中值滤波以去除单独的刺和洞；然后移除10％的边界点，这些边界点到质心的距离与平均距离相差最大；最后用立方的b样条拟合剩下的边界点，形成精分割。
[0109] 以下对本发明的实施例的实验与结果作进一步详细说明。
[0110] 在基于ITK(Insight Segmentation and Registration Toolkit，是美国国家卫生院下属的国立医学图书馆开发的一款医学
图像处理软件包)的开源配准工具包elastix上实现加权α-MI度量。所有的配准实验都是在64位Windows系统(Intel 48-cores 2.60GHz CPU和128GB内存)上执行。DCNN模型的训练实验是在64位Windows系统(Intel 8-cores 3.40GHz CPU,NVIDIA GeForce GTX 1080和16GB内存)上执行。
[0111] (一)数据集和验证方法
[0112] 实验采用17例患者的T2-加权MRI宫颈图像数据，在五周的时间里，每位患者每周接受一次扫描。图像的大小为512×512×30，体素大小为0.625×0.625×4.5毫米。在配准之前，所有图像都被裁剪为210×250×30大小。由两名临床专家对宫颈(cervix)、膀胱(bladder)和直肠(rectum)进行手动分割，它们将在
算法性能评估中被视为金标准。
[0113] 将数据分为训练集和测试集。训练集包括所有患者的第一周和第二周图像，用于调整配准参数和训练DCNN模型。剩下的51幅图像组成测试集用于方法的验证。在atlas-to-target配准实验中，将aMI-MST和aMI-WMST与LMI(Localized Mutual Information,局部互信息)进行比较，测量各算法自动分割与金标准之间的DSC(Dice Similarity Coefficient,Dice相似系数)。双侧Wilcoxon试验中p＜0.05的值被认为是有统计学意义的差异。
[0114] (二)参数设置
[0115] 对训练集通过试错的方式选择配准参数，在所有可形变配准之前还需执行基于MI的刚性配准。针对可形变配准，一种基于高斯平滑的多分辨率策略被使用。在x和y方向使用σ＝4.0，2.0，1.0的范围，在z方向使用σ＝2.0，1.0，0.5的范围。三种分辨率分别采用网格间距为80mm、40mm和20mm的参数化b样条变换模型。
[0116] 设置A＝50,τ＝0.602,a＝2000,和600次迭代以优化aMI-MST，随机样本的数量设置为N＝5000。对于aMI-WMST和LMI，在每次迭代中随机选取大小为403立方毫米的局部邻域，设置A＝200,τ＝0.6,a＝2000,和2000次迭代，随机样本的数量设置为N＝2000。
[0117] 在MST的实现过程中，选择一个bucket大小为50、误差边界值为10.0的kD树。其中，k＝20。此外，α值被设置为0.99。在VGG-19模型的迁移学习中，学习率设置为0.0001，批量大小设置为20。
[0118] (三)粗分割结果分析
[0119] 为了对测试集的34幅图像进行粗分割，进行了同一患者不同周次图像之间的配准(共17×2次)。对每位患者，选择某一周的图像作为模板图像，而下周的图像作为目标图像。通过主成分分析，将多维笛卡尔特征集压缩成六个分量的子集合。三种方法的Dice系数箱线图如图4所示。结果表明，LMI的分割精度明显低于其他方法。与aMI-MST相比，aMI-WMST的-3
中值重叠率明显增加，其中宫颈部位由0.834增长到0.843(p＝1.0×10 )，膀胱部位由
0.804增长到0.826(p＝1.2×10-2)，直肠部位由0.767增长到0.782(p＝1.1×10-3)。
[0120] 如图4所示，不同方法在不同解剖结构下的DSC箱线图。对于每个结构，星号表示中值重叠率与前一列相比存在统计显著性差异。
[0121] (四)精分割结果分析
[0122] 在使用aMI-WMST进行可形变配准之后，34幅目标图像的宫颈和膀胱结构被微调的VGG-19模型进行进一步的提炼。与粗分割相比，精分割的中值重叠率显著增加，其中宫颈部位由0.843增加至0.882，膀胱部位由0.826增加至0.857。图5为分割结果对比展示图，其中，(a)为金标准的分割；(b)为粗分割结果；(c)为最终分割结果。虽然在最终结果中可以发现膀胱和宫颈附近软组织的复杂性，但是通过精分割得到的边界比粗分割更接近于金标准。
[0123] 显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些
修改和变型属于本发明
权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
