CN109215064B - 一种基于超像素向导的医学图像配准方法 - Google Patents

Abstract

本发明公开了一种基于超像素向导的医学图像配准方法，包括步骤：1)对参考图像和浮动图像采用基于B样条的FFD配准算法进行刚性预配准，获得预配准后的参考图像和浮动图像；2)对预配准后的参考图像和浮动图像，进行下采样，形成多分辨率金字塔，级数为l；3)多分辨率方案的图像配准从最粗的分辨率开始，由粗至精的逐步配准，在当前分辨率l下，经迭代后，对已获得的空间变换u_l进行上采样得到u_l‑1，然后利用空间变换u_l‑1对金字塔第l‑1级的图像进行变换，将获得的图像再次进入迭代，直到金字塔图像最低层结束；4)将最终获得的空间变换u作用于刚性配准后的浮动图像，获得最终的配准后结果图像。本发明获得了优于基于B样条的FFD模型、DD模型等当前主流DIR配准模型的配准精度和鲁棒性。

Description

一种基于超像素向导的医学图像配准方法

技术领域

本发明涉及医学图像处理及应用的技术领域，尤其是指一种基于超像素向导的医学图像配准方法。

背景技术

肺癌是当今全世界发病率和死亡率增长最快、对人类健康和生命威胁最大的恶性肿瘤之一。肺癌切除术后胸部随访CT(Computed Tomography,CT)图像之间的配准是预测术后肺功能、监测术后肿瘤的复发和转移情况、评估术后残肺代偿效果以及判断预后并指导治疗等的关键指标。然而，由于肺部固有的易变形性，再加上器官蠕动、呼吸运动及患者移动等特征，使得胸部随访图像之间存在复杂大形变，对它们精确地配准仍然是一项艰巨的任务。近年来，针对肺部的形变图像配准(Deformation Image Registration，DIR)已经取得了重要进展：该类算法比以往更快，更精确，更易于使用。然而，DIR通常依赖于所寻求的变换是平滑的这一假设，从而阻止对那些不连续性解剖变化(例如组织生长/收缩或肿瘤反应)的正确估计。因此，为了能够处理由肺部代偿生长、肿瘤组织收缩/膨胀以及肿瘤反应等原因引起的复杂大形变，对肺部的生理特性、运动特点和解剖结构等先验知识来向导DIR模型优化进行研究具有重要的理论意义和临床应用价值。

当前，Demons模型是肺部图像进行形变配准最流行使用的一类DIR模型。然而，该类DIR模型假设所寻求的变换是平滑的，从而阻止对那些不连续性解剖变化(例如肺部膨胀/收缩、体重减轻/增加或肿瘤反应)的正确估计。因此，为了能够处理由肺部代偿生长、肿瘤组织收缩/膨胀以及肿瘤反应等原因引起的复杂大形变，需要结合肺部的生理特性、运动特点和解剖变化特点等先验知识来向导DIR模型优化，从而提高肺部图像配准的精度和鲁棒性。近年来，一些学者对向导DIR模型进行了研究。例如，Kim等人将人工勾画的器官轮廓和一些特殊点作为约束项结合于基于B样条的FFD模型，目的是建立临床患者数据的基准配准。Kearney等人提出了一种基于标记向导的图像配准(Landmark-Guided ImageRegistration，LGIR)算法，该算法采用梯度相似性度量来建立标记以帮助向导Demons模型的目标函数趋向收敛。Wu等人提出了一种基于生长轨迹或时空对应关系向导的形变配准算法，该算法解决了不同时期婴儿大脑图像的配准难题，并且与现有最先进的配准算法相比产生了更高的配准精度。Gu等人提出了一种基于轮廓向导的形变配准算法(Contour-Guided Deformable Image Registration，CG-DIR)，以提高自适应放射治疗的变形向量场的准确性和一致性。该算法使用勾画的结构集对之间的灰度匹配项对原始Demons模型中的目标函数进行正则化，来实现将医生勾画的轮廓线结合进配准模型。然而，不同医生可能产生不同的勾画结果，甚至同一医生不同时间点都可能勾画出不同的轮廓，因此配准结果的精度和一致性受医生主观经验影响很大，并且人工勾画过程耗时较长，对大数据集进行分析可能是非常费力的。

本发明提供一种基于超像素向导的医学图像配准方法，将肺部的生理特性、运动特点和解剖结构等先验知识结合进微分同胚Demons(Diffeomorphic Demons，DD)模型的目标函数中，构建了一种新配准算法模型。新模型将待配准两图像的超像素二值表示之间距离作为新匹配项以向导基于DD模型的目标函数趋向收敛，克服了DD算法仅依赖灰度信息驱动形变场引起的形变不足问题，有助于提高由边界运动或解剖变化所引起复杂形变的配准精度，同时增强配准鲁棒性。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提出了一种基于超像素向导的医学图像配准方法，获得了优于基于B样条的FFD模型、DD模型等当前主流DIR配准模型的配准精度和鲁棒性，对大形变区域和解剖变化区域的精确匹配具有更重要的临床应用意义和参考价值。

为实现上述目的，本发明所提供的技术方案为：一种基于超像素向导的医学图像配准方法，包括以下步骤：

1)对参考图像和浮动图像采用基于B样条的FFD配准算法进行刚性预配准，获得预配准后的参考图像和浮动图像；

2)对预配准后的参考图像和浮动图像，进行下采样，形成多分辨率金字塔，级数为l；

3)多分辨率方案的图像配准从最粗的分辨率开始，由粗至精的逐步配准，在当前分辨率l下，经如下迭代过程之后，对已获得的空间变换u_l进行上采样得到u_l-1，然后利用空间变换u_l-1对金字塔第l-1级的图像进行变换，将获得的图像再次进入如下迭代过程，直到金字塔图像最低层结束；其中，所述迭代过程如下：

3.1)分别计算当前分辨率l下参考图像和浮动图像的超像素的二值表示B₀(x)和B₁(x)；

3.2)将当前超像素二值表示B₀(x)和B₁(x)之间的欧氏距离构造的一个新的匹配项结合于经典DD模型的目标能量函数中，构建新目标能量函数；

3.3)根据新目标能量函数，利用空间变换u_l-1与更新速度场v_l-1关系式和新目标能量函数的一阶最优性条件，获取使得新目标函数最小化的空间变换函数；

3.4)判断是否达到最大迭代次数，若没有达到，则返回步骤3.1)，继续迭代；否则判断是否达到最大分解级数，若没有达到，则对低分辨率下得到的形变矩阵进行超采样，将其作为上一级高分辨率的初始变换，返回步骤3.1)，否则配准结束，进入步骤4)；

4)将最终获得的空间变换u作用于刚性配准后的浮动图像，获得最终的配准后结果图像。

在步骤3.1)中，所述参考图像和浮动图像的超像素的二值表示计算，包括以下步骤：

3.1.1)初始化聚类种子点，对于一幅有N个像素点的图像I进行SLIC超像素分割，分割块数根据分割对象人为设定，这也是SLIC算法的突出优势之一；如果设定分割的超像素个数为K，即整幅图像被均匀分割为K个网格，那么每一网格内含有N/K像素；网格中心即为初始化的聚类种子点，则相邻超像素之间的步长近似为

3.1.2)利用SLIC算法分割获取图像超像素的标签图B(x)；

3.1.3)遍历超像素标签图中每一标签，基于梯度边缘检测算法判断是否属于超像素边缘，如果是则将其标记为255，即B(x)＝255，否则标记为0，即B(x)＝0，获得了超像素的二值表示B(x)。

在步骤3.2)中，在DD算法模型基础上，采用新目标能量函数：

式中，V₁和V₀分别表示参考图像和浮动图像的灰度值，B₁和B₀分别表示参考图像和浮动图像的超像素二值表示；γ和β分别表示新添加匹配项和变换函数光滑程度在目标能量函数中所占权重；u为空间变换；c作为辅助的空间变换，目的是简化优化过程，c和u存在关系式

其中v表示更新速度向量场；第三项为形变场的偏差项；

在步骤3.3)中，对新目标能量函数，利用其一阶最优性条件，获取使得新目标函数最小化的空间变换函数过程中，包括以下步骤：

3.3.1)空间变换u和更新速度场v的关系式为

3.3.2)将E(u)的三项分别进行一阶泰勒展开，得到如下近似表达式：

exp(v)-Id≈v

式中，J_V ^P和J_B ^P分别表示第一项和第二项的空间雅克比矩阵；

3.3.3)利用E(u)一阶最优性条件

获得方程的解，即为更新速度场，即为更新速度场：

3.3.4)采用关系式u＝exp(v)得到最优的空间变换函数u。

本发明与现有技术相比，具有如下优点与有益效果：

1、本发明首先采用基于像素灰度的刚性配准算法对图像进行全局粗配准，使得参考图像和浮动图像达到空间位置上一致，然后以粗配准结果作为初值，再采用所提SG-DIR算法在多分辨率策略的框架中进行精配准，提高了配准精度，降低了计算复杂度。

2、本发明采用待配准两图像的超像素二值表示之间的距离作为新匹配项以向导DD模型中的目标函数趋向收敛，不仅在小形变区域获得了高精度的配准结果，而且在大形变区域和解剖变化区域均获得了高精度的配准结果，在临床中的实际应用有着重要的参考价值。

3、本发明相比于FFD和DD模型的形变配准算法，获得了更快的收敛速率、更高的配准精度和配准效率。

附图说明

图1为本发明逻辑流程示意图。

图2为本发明由一幅CT序列片图像肺实质分割前后的超像素二值表示图。

图3为本发明不同配准算法对小形变图像之间的配准结果分析。

图4为本发明不同配准算法对大形变图像之间的配准结果分析。

图5为本发明不同配准算法对存在解剖变化的图像之间配准结果分析。

图6为本发明客观度量均值随迭代次数变化的曲线图。

具体实施方式

下面结合具体实施例对本发明作进一步说明。

如图1至图6所示，本实施例所提供的基于超像素向导的医学图像配准方法，使用的实验平台是Matlab7.11，计算机主机配置：Intel(R)处理器，CPU主频3.0GHz，内存4.0GB，其包括以下步骤：

1)对参考图像和浮动图像采用基于B样条的FFD配准算法进行刚性预配准，获得预配准后的参考图像和浮动图像，确保在形变配准过程中的参考图像和浮动图像是空间对齐的。

2)对预配准后的参考图像和浮动图像，进行下采样，形成多分辨率金字塔，级数为l，其中l＝1,2,3。

3)多分辨率方案的图像配准从最粗的分辨率开始，由粗至精的逐步配准。在当前分辨率l下，经如下迭代过程之后，对已获得的空间变换u_l进行上采样得到u_l-1。

多分辨率方案中，从高分辨率到低分辨率的最大迭代次数分别为100次、50次和20次，目标能量函数前后两次迭代的误差值不大于1×10^-4为迭代停止准则。多分率策略避免了局部极值的产生,有助于提高配准鲁棒性,并降低计算代价。

利用空间变换u_l-1对金字塔第l-1级的图像进行变换，将获得的图像再次进入如下迭代过程，直到金字塔图像最低层结束。

3.1)计算当前分辨率l下图像的超像素二值表示B(x)，包括以下步骤：

3.1.1)初始化聚类种子点。对于一幅有N个像素点的图像I进行SLIC超像素分割，分割块数可根据分割对象人为设定，这也是SLIC算法的突出优势之一。如果设定分割的超像素个数为K，即整幅图像被均匀分割为K个网格，那么每一网格内含有N/K像素。网格中心即为初始化的聚类种子点，则相邻超像素之间的步长近似为

3.1.2)利用SLIC算法分割获取图像超像素的标签图B(x)。

3.1.3)遍历超像素标签图中每一标签，基于梯度边缘检测算法判断是否属于超像素边缘，如果是则将其标记为255，即B(x)＝255，否则标记为0，即B(x)＝0，获得了超像素的二值表示B(x)。

超像素的二值表示能够精确反应和描述边界运动和解剖变化，将其作为先验信息来向导具有复杂形变图像的配准过程，有助于改善配准精度、加速能量函数收敛以及提高配准的鲁棒性。

3.2)将当前参考图像和浮动图像的超像素二值表示B₁(x)和B₀(x)之间的欧氏距离构造的匹配项结合于经典DD模型的目标能量函数中，构建新目标能量函数：

式中，V₁和V₀分别表示参考图像和浮动图像的灰度值，B₁和B₀分别表示参考图像和浮动图像的超像素二值表示；γ和β分别表示新添加匹配项和变换函数光滑程度在目标能量函数中所占权重；u为空间变换；c作为辅助的空间变换，目的是简化优化过程，c和u存在关系式

其中v表示更新速度向量场；第三项为形变场的偏差项。

新目标能量函数结合了一个由超像素二值表示之间的欧氏距离构造的新匹配项

将其作为先验信息来向导具有复杂形变图像的配准过程，有助于改善配准精度、加速能量函数收敛以及提高配准的鲁棒性。

3.3)根据新目标能量函数，利用空间变换u_l-1与更新速度场v_l-1关系式和新目标能量函数的一阶最优性条件，获取使得新目标函数最小化的空间变换函数的过程中，包括以下步骤：

3.3.1)空间变换u和更新速度场v的关系式为

3.3.2)将E(u)的三项分别进行一阶泰勒展开，得到如下近似表达式：

exp(v)-Id≈v

式中，J_V ^P和J_B ^P分别表示第一项和第二项的空间雅克比矩阵。

3.3.3)利用E(u)一阶最优性条件

获得方程的解，即为更新速度场，即为更新速度场：

3.3.4)采用关系式u＝exp(v)得到最优的空间变换函数u。

所提出的更新速度场v避免了仅依靠灰度差和梯度信息的弊端，增加了边界运动及解剖变化的先验信息的平滑约束，有助于加速收敛速度，减少形变配准误差，增强配准的鲁棒性。

为了分析配准性能，本实施例以11例术后随访检查的胸部CT为配准实验对象，采用DSC度量将所提SG-DIR算法与经典DD算法进行比较，获得的DSC度量值如下表1所示：

表1所提算法和经典DD算法的DSC度量值

从表1中可以看出，所提SG-DIR算法获得的DSC度量值均高于经典DD算法，这表明所提SG-DIR算法的配准性能优于经典DD算法，这与主观可视化图像的观察结果一致，因而所提算法对大形变和灰度变化的配准性能较好。

为了进一步综合的评估所提算法的配准性能，本实施例对所有已给定的11组待配准图像对分别采用FFD算法、经典DD算法和所提算法进行配准，输出每组配准结果的客观定量度量值及配准完成时间；然后，对它们取均值，分别由如下符号表示：mRSSD、mDSC、mNCC、mRMSE和mT/s；最后，将这些客观度量及时间均值显示在同一表格内，以便比较分析，如下表2所示。

表2不同配准算法的客观评估度量值

从表2中可以看出，所提SG-DIR算法获得的mRSSD和mRMSE的值最小，同时获得的mDSC和mNCC值最大，配准时间与DD算法相近，所有客观度量均值表明所提算法在灰度匹配和形变程度之间达到最好平衡，配准误差最小，配准后结果图像与参考图像之间重叠度和相关性最优，并且也保持了较高的配准效率。

3.4)判断是否达到最大迭代次数，若没有达到，则返回步骤3.1)，继续迭代；否则判断是否达到最大分解级数，若没有达到，则对低分辨率下得到的形变矩阵进行超采样，将其作为上一级高分辨率的初始变换，返回步骤3.1)，否则配准结束，进入步骤4)；

4)将最终获得的空间变换u作用于刚性配准后的浮动图像，获得最终的配准后结果图像。

综上所述，在采用以上方案后，本发明增加了基于边界运动及解剖变化先验信息的新匹配项，避免了仅依靠灰度差和梯度信息的弊端，减少了陷入局部极值的风险，改善了配准精度，增强了配准鲁棒性，对在临床中的实际应用有着重要的参考价值。

以上所述实施例只为本发明之较佳实施例，并非以此限制本发明的实施范围，故凡依本发明之形状、原理所作的变化，均应涵盖在本发明的保护范围内。

Claims (2)

1.一种基于超像素向导的医学图像配准方法，其特征在于，包括以下步骤：

1)对参考图像和浮动图像采用基于B样条的FFD配准算法进行刚性预配准，获得预配准后的参考图像和浮动图像；

2)对预配准后的参考图像和浮动图像，进行下采样，形成多分辨率金字塔，级数为l；

3)多分辨率方案的图像配准从最粗的分辨率开始，由粗至精的逐步配准，在当前分辨率l下，经如下迭代过程之后，对已获得的空间变换u_l进行上采样得到u_l-1，然后利用空间变换u_l-1对金字塔第l-1级的图像进行变换，将获得的图像再次进入如下迭代过程，直到金字塔图像最低层结束；其中，所述迭代过程如下：

3.1)分别计算当前分辨率l下参考图像和浮动图像的超像素的二值表示B₀(x)和B₁(x)；

3.2)将当前超像素二值表示B₀(x)和B₁(x)之间的欧氏距离构造的一个新的匹配项结合于经典DD模型的目标能量函数中，构建新目标能量函数；

在DD算法模型基础上，采用新目标能量函数：

式中，V₁和V₀分别表示参考图像和浮动图像的灰度值，B₁和B₀分别表示参考图像和浮动图像的超像素二值表示；γ和β分别表示新添加匹配项和变换函数光滑程度在目标能量函数中所占权重；u为空间变换；c作为辅助的空间变换，目的是简化优化过程，c和u存在关系式

其中v表示更新速度向量场；第三项为形变场的偏差项；

3.3)根据新目标能量函数，利用空间变换u_l-1与更新速度场v_l-1关系式和新目标能量函数的一阶最优性条件，获取使得新目标函数最小化的空间变换函数；

对新目标能量函数，利用其一阶最优性条件，获取使得新目标函数最小化的空间变换函数过程中，包括以下步骤：

3.3.1)空间变换u和更新速度场v的关系式为

3.3.2)将E(u)的三项分别进行一阶泰勒展开，得到如下近似表达式：

exp(v)-Id≈v

式中，J_V ^P和J_B ^P分别表示第一项和第二项的空间雅克比矩阵；

3.3.3)利用E(u)一阶最优性条件

获得方程的解，即为更新速度场，即为更新速度场：

3.3.4)采用关系式u＝exp(v)得到最优的空间变换函数u；

3.4)判断是否达到最大迭代次数，若没有达到，则返回步骤3.1)，继续迭代；否则判断是否达到最大分解级数，若没有达到，则对低分辨率下得到的形变矩阵进行超采样，将其作为上一级高分辨率的初始变换，返回步骤3.1)，否则配准结束，进入步骤4)；

4)将最终获得的空间变换u作用于刚性配准后的浮动图像，获得最终的配准后结果图像。

2.根据权利要求1所述的一种基于超像素向导的医学图像配准方法，其特征在于：在步骤3.1)中，所述参考图像和浮动图像的超像素的二值表示计算，包括以下步骤：

3.1.1)初始化聚类种子点，对于一幅有N个像素点的图像I进行SLIC超像素分割，分割块数根据分割对象人为设定，这也是SLIC算法的突出优势之一；如果设定分割的超像素个数为K，即整幅图像被均匀分割为K个网格，那么每一网格内含有N/K像素；网格中心即为初始化的聚类种子点，则相邻超像素之间的步长近似为

3.1.2)利用SLIC算法分割获取图像超像素的标签图B(x)；

3.1.3)遍历超像素标签图中每一标签，基于梯度边缘检测算法判断是否属于超像素边缘，如果是则将其标记为255，即B(x)＝255，否则标记为0，即B(x)＝0，获得了超像素的二值表示B(x)。

Images