CN108364307A - 一种基于α-Renyi互信息的医学图像配准方法 - Google Patents

Abstract

本发明公开了一种基于α‑Renyi互信息的医学图像配准方法(Medical Image Registration based on the α‑Renyi mutual information，RMIR)，包括如下步骤：获取图像像素信息步骤，获取参考图像和浮动图像的像素信息；获取图像灰度信息步骤，获取所述参考图像和所述浮动图像的灰度信息；获取图像互信息量步骤，获取所述参考图像和所述浮动图像的Renyi互信息量，图像配准步骤，将所述获取图像互信息量步骤中获取的所述Renyi互信息量进行刚性配准，通过所述配准函数获得最佳变换，使得所述Renyi互信息达到最大值。通过Renyi熵引入互信息量的度量中，并通过配准函数进行配准，有效提高了配准精度，并通过大量配准实验，找寻合适的α初始参数，有效降低了配准过程中的运算负载。

Description

一种基于α-Renyi互信息的医学图像配准方法

技术领域

本发明涉及一种图像处理方法，尤其涉及一种医学图像配准方法。

背景技术

在为患者诊断过程中，经常要将同一患者多幅图像放在一起分析，得到 该患者的多方面的综合信息，从而提高医学诊断准确率，于是需要将多幅不 同时间、不同场景、不同模态的图像进行空间几何变换，使得代表相同结果 的像素(体素)在空间上达到一致，通常有基于图像特征和基于图像灰度的 配准方法，其中基于图像灰度的配准方法利用图像灰度信息，配准精度较高， 例如Shannon互信息进行配准(Medical Image Registrationbased on the Shannonmutual information，SMIR)，但是由于互相信息函数包含了大量的局 部极小值，容易陷入局部最优，会不断循环迭代地搜索解空间，计算量巨大， 并且如果初始值选择不当，会消耗更多的搜索时间，甚至造成配准失败。

发明内容

本发明的目的是提供一种基于α-Renyi互信息的医学图像配准方法，通 过Renyi熵引入互信息量的度量中，并通过配准函数进行配准，有效提高了 配准精度，并通过大量配准实验，找寻合适的初始参数，有效降低了配准过 程中的运算负载。

本发明公开了一种基于α-Renyi互信息的医学图像配准方法，包括如下 步骤：获取图像特征步骤，获取参考图像和浮动图像的图像特征；获取图像 灰度信息步骤，获取所述参考图像和所述浮动图像的灰度信息；获取图像互 信息量步骤，通过所述参考图像和所述浮动图像的所述图像特征和灰度信息 获取Renyi互信息量，所述Renyi互信息量由下式获得：

RMI(R,F)＝H_α(R)+H_α(F)-H_α(R,F)

其中，RMI(R，F)为参考图像R和浮动图像F的Renyi互信息量，H_α(R) 和H_α(F)分别为R和F的Renyi熵，H_α(R,F)是图像R和F的Renyi熵；所 述Renyi熵由下式获得：

图像配准步骤，将所述获取图像互信息量步骤中获取的所述互信息量进 行刚性配准，配准函数如下式：

其中，T(F)表示对浮动图像进行刚性变换；

通过所述配准函数获得最佳变换，使得所述Renyi互信息达到最大值。

进一步地，上述的基于α-Renyi互信息的医学图像配准方法中所述Renyi 熵中参数α的取值范围为：-0.1<α<0。

进一步地，上述的基于α-Renyi互信息的医学图像配准方法中所述Renyi 熵中参数α的取值范围为：0<α<1。

进一步地，上述的基于α-Renyi互信息的医学图像配准方法中所述Renyi 熵中参数α的取值范围为：1<α<1.3。

通过引入Renyi熵进行互信息量计算，避免了在配准过程运算负载过重 等问题，提高了配准精度。

附图说明

图1是本实施例所示方法的流程图；

图2a-c是本实施例的参考图像；

图3是α<0时，Renyi互信息RMI变化曲线图；

图4是α<0时，Renyi互信息运算负荷示意图；

图5a-b是0<α<1时，Renyi互信息RMI变化曲线图；

图6是0<α<1时，Renyi互信息运算负荷示意图；

图7是α>1时，Renyi互信息RMI变化曲线图；

图8是α>1时，Renyi互信息运算负荷示意图；

图9a-c是α＝1.1和α＝0.95时，对图2a-c进行RMI和SMI配准结果变 化示意图；

图10是α＝1.1和α＝0.95时，对图2a-c进行RMI和SMI配准运算符合 示意图；

图11a-b为单模态下对图2a-c变化后的CT浮动图像；

图12a-b为单模态下对图2a-c变化后的MR浮动图像；

图13a-b为单模态下对图2a-c变化后的PET浮动图像；

图14a-b为实验1的CT图像配准结果示意图；

图15a-b为实验1的MR图像配准结果示意图；

图16a-b为实验1的PET图像配准结果示意图；

图17a-b为实验2的CT图像配准结果示意图；

图18a-b为实验2的MR图像配准结果示意图；

图19a-b为实验2的PET图像配准结果示意图；

图20a-b为实验3的CT图像配准结果示意图；

图21a-b为实验3的MR图像配准结果示意图；

图22a-b为实验3的PET图像配准结果示意图；

图23a-b为多模态下对图2a-c变化后的CT浮动图像；

图24a-b为多模态下对图2a-c变化后的MR浮动图像；

图25a-b为多模态下对图2a-c变化后的PET浮动图像；

图26a-f为实验4的图像配准结果示意图；

图27a-f为实验5的图像配准结果示意图；

图28a-f为实验6的图像配准结果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施 方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例 性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结 构和技术的描述，以避免不必要地混淆本发明的概念。

本实施例所示的一种基于α-Renyi互信息的医学图像配准方法(Medical ImageRegistrationbased on theα-Renyi mutual information，缩写为：RMIR)， 包括如下步骤：获取图像像素信息步骤、获取图像灰度信息步骤、获取图像 互信息量步骤和图像配准步骤，

获取图像特征步骤，获取参考图像和浮动图像的图像特征；

获取图像灰度信息步骤，获取所述参考图像和所述浮动图像的灰度信息；

获取图像互信息量步骤，通过所述参考图像和所述浮动图像的所述图像 特征和灰度信息获取Renyi互信息量，所述Renyi互信息量由下式获得：

RMI(R,F)＝H_α(R)+H_α(F)-H_α(R,F)

其中，RMI(R，F)为参考图像R和浮动图像F的Renyi互信息量，H_α(R) 和H_α(F)分别为R和F的Renyi熵，H_α(R,F)是图像R和F的Renyi联合熵； 所述Renyi熵由下式获得：

图像配准步骤，将所述获取图像互信息量步骤中获取的所述互信息量进 行刚性配准，配准函数如下式：

通过所述配准函数获得最佳变换，使得所述Renyi互信息达到最大值。

本实施例中基于α-Renyi互信息的医学图像配准方法中所述Renyi熵中 参数α的取值范围为：-0.1<α<0。

本实施例中基于α-Renyi互信息的医学图像配准方法中所述Renyi熵中 参数α的取值范围为：0<α<1。

本实施例中基于α-Renyi互信息的医学图像配准方法中所述Renyi熵中 参数α的取值范围为：1<α<1.3。

下面对本实施例的工作原理进行说明。

本实施例中，参考图像R和浮动图像F大小均调整为M×N，L级灰度， 左上角像素点位置为(1,1)，参考图像R中点(x，y)的灰度值为r(x，y)， 对应的浮动图像F中点(x，y)的灰度值为f(x，y)，于是可以用图像的Shannon 熵来度量图像中灰度分布的聚集特征所包含的信息量，从而可以得出：

其中，p(k)为灰度值k出现的概率。

医学图像配准过程中，当参考图像R和浮动图像F的空间位置一致时， 如对应像素的灰度互信息达到最大值，即一幅图像表达的关于另外一幅图像 的信息最多，可以将这个称为互信息，可以通过下式计算：

SMI(R,F)＝H(R)+H(F)-H(R,F) (2)

其中，SMI(R，F)为参考图像R和浮动图像F的Shannon互信息量， H(R)和H(F)分别为参考图像R和浮动图像F的熵，H(R，F)是参考 图像R和浮动图像F的联合熵，则基于互信息量的刚性配准可以表示为：

于是需要寻找一个最佳变换T_optimal，使得互信息量达到最大。在二维空间， 医学图像刚体变换实质上包括了平移和旋转变换，即：

其中，Δx，Δy，φ分别为垂直平移量、水平平移量和旋转角度。

由于在对式(4)求解最优解时，由于互信息量函数包含了大量的局部极 小值，非常容易陷入局部最优，于是在求解过程中，会不断循环迭代搜索解 空间，计算量非常大，同时如果初始值选择不当，会消耗更多搜索空间，甚 至造成配准失败。

Shannon熵为α阶广义熵，其中α阶广义熵的定义如下式：

其一阶广义熵即为Shannon熵，由此出发，可以得出图像的Renyi熵， 于是可以基于Renyi互信息量的图像配准可以解决该问题，首先定义Renyi 熵如下所示：

其中，α(α≠1)是一个正实数，p(k)为灰度值k出现的概率，代 入公式(2)，用Renyi熵代替Shannon熵，则公式(2)变为：

RMI(R,F)＝H_α(R)+H_α(F)-H_α(R,F) (7)

其中，RMI(R，F)为参考图像R和浮动图像F的Renyi互信息量，H_α(R) 和H_α(F)分别为R和F的Renyi熵，H_α(R,F)是图像R和F的Renyi联合熵。

RMI为R和F的配准函数，则基于Renyi互信息量的刚性配准可以表示 为：

其中，T(F)表示对浮动图像进行刚性变换；

通过公式(7)便可以得出一个最佳变换T_optimal，使得Renyi互信息量达 到最大，从而完成医学图像配准。

下面通过美国Vanderbilt大学RREP(Retrospective Registration EvaluationProjection)项目组的国际通用刚性配准脑部图像数据，CPU为 Dual-CoreE5500+2.5GHz，内存为2GB，操作系统为Windows XP并采用 Matlab7.1为实验软件进行医学图像配准，并采用单纯形法作为多参数寻优方 法。

参照图2a、图2b、图2c，分别为：CT(Computed Tomography，即电 子计算机断层扫描)第1幅图像、MR T1(Magnetic Resonance，即磁共振 成像)第1幅图像和PET(PositronEmission Computed Tomography，即正 电子发射型计算机断层显像)第2幅图像作为参考图像，大小分别为512× 512、256×256和128×128。

为了获得对比，下面对α-Renyi互信息和Shannon互信息均进行了实验。 本实施例中，参考图像R和浮动图像F大小均调整为M×N，L级灰度。

其中，一次Shannon互信息量运算负荷如下表所示：

表1一次Shannon互信息量运算负载

其中，一次α-Renyi互信息量运算负荷如下表所示：

表2一次α-Renyi互信息量运算负载

从表1和表2中可知，一次Shannon互信息量运算负载和一次α-Renyi 互信息量运算负载分别为：(3MN+4L2+8L)、(3MN+2L2+ L2α²+2Lα+4L+3)，由于图像大小、图像灰度级为确定值，因此参数α的合 适选择对α-Renyi互信息量运算负载有重大影响。为了进一步讨论参数α的 选取范围，选择CT图像作为参考图像(图2a所示)，并以α-Renyi互信息 量作为相似性测度。由于图像角度旋转在图像配准时非常耗时，而且精度要 求高，因此，我们通过图像角度旋转对参考图像进行自配准来选择合适的参 数α。

1.α<0

参考图3，在实验中，通过图像角度旋转参考图像进行自配准，参数 α∈{-0.02,-0.04,-0.05,-0.06,-0.08,-0.1,-0.2,-0.3,-0.4,-0.5,-0.6,-0.7,-0.8,-0.9}，取不同值时RMI(公式(7))变化曲线分别如图2所示，可见当-0.1<α<0时，RMI曲 线较为光滑，能得到精确的配准参数，并且不同的α参数对应的RMI曲线差 异很小，曲线几乎重合；当α<-0.1时，RMI曲线出现了非“零点”的局部极 值点，α越小，非“零点”极值点越多，显然不利于获得正确的配准参数， 易陷入局部极值点，因此α<-0.1不适合用于配准。参考图4，α＝-0.1时，配 准运算负载最低，其它α参数配准运算负载相对较高。

2.0<α<1

在实验中，通过图像角度旋转参考图像进行自配准，参数α分别取 α∈{0.02,0.04,0.05,0.06,0.08,0.1,0.2,0.3,0.4,0.5}和 α∈{0.6,0.7,0.8,0.9,0.92,0.94,0.96,0.98}这两个集合，取不同值时RMI变化曲线 分别如图5a和图5b所示，可见当0<α<1时，RMI曲线较为光滑，均能得到精 确的配准参数，α越大，RMI曲线越平滑，越利于获得正确的配准参数；当 0<α<0.1和0.9<α<1.0时，RMI曲线形状分别与α＝0.1和α＝0.9相似且几乎重合，也能得到精确的配准参数，如图6所示。整体上看，配准运算负载差异 不大。

3.α>1

参考图7，在实验中，通过图像角度旋转参考图像进行自配准，参数α分 别取α∈{1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2,2.1,2.2,2.3,2.4,2.5,3}，取不同值时RMI变化曲线分别如图7所示，由于上述参数得到的RMI曲线差异很小，无法显 示全部RMI曲线，因此图7只显示了部分α对应的RMI曲线，此时均能得 到精确的配准参数，同时当α>1能得到精确的配准参数，而当1<α<1.3时， RMI曲线较为平滑，α越来越大时，RMI曲线在零点附近越陡峭，曲线变得 不平滑；当α＝2时，运算负载最低。

可以见得，当-0.1<α<0,0<α<1和1<α<1.3均能用于图像配准，且RMI 曲线平滑，从而可以降低运算负载。接着选取α＝1.1和α＝0.95作为RMI的 α参数值，对图2a、图2b和图2c分别使用RMI和SMI进行配准，图9a、 图9b、图9c和图10为实验结果。可以看出，当α＝1.1和α＝0.95时，RMI 和SMI形状很相似，配准运算负载也非常接近。

接着将图2a、图2b和图2c作为参考图像，对每幅参考图像做如表3所 示的变换，得到浮动图像，分别如图11a和图11b、图12a和图12b、图13a 和图13b所示，

表3单模态浮动图像变换参数

为了比较配准精度，定义误差η_i为：

其中，Δi_s表示浮动图像相对于参考图像的标准变换参数；Δi表示通过 配准获取的浮动图像变换参数，于是总误差η为：

先进行单模态图像模拟试验，在α＝0.95时，将RMIR(Medical ImageRegistration based on theα-Renyi mutual information，即基于α-Renyi互 信息的医学图像配准方法)，用于单模态配准得出实验结果1，参考图14a和 图14b、图15a和图15b、图16a和图16b，以及表4，

表4实验1单模态配准性能

在α＝1.1时，将RMIR用于单模态配准得出实验结果2，参考图17a和 图17b、图18a和图18b、图19a和图19b，以及表5，

表5实验2单模态配准性能

将SMIR(Medical Image Registration based on the Shannon mutualinformation，即基于Shannon互信息的医学图像配准方法)用于单模态配准 得出实验结果3，参考图20a和图20b、图21a和图21b、图22a和图22b， 以及表6，

表6实验3基于互信息量单模态配准性能

从以上实验1、实验2和实验3结果来看，RMIR和SMIR相比较，在6个配 准实验中，在配准速度方面，SMIR、RMIR(0.95)和RMIR(1.1)三者非常接近， 差距均在0.5秒内。因此，从整体上来看，RMIR配准速度和SMIR相当。

从以上实验1、实验2和实验3结果来看，RMIR和SMIR相比较，在CT+CT1 配准时，SMIR配准速度比RMIR快，但RMIR接近SMIR配准速度；在CT+CT2 配准时，RMIR(1.1)配准速度有较强优势，RMIR(0.95)耗时仍然最长；在 MR+MR1配准时，RMIR(1.1)配准最快，RMIR(1.1)耗时最长；在MR+MR2 配准稍时，RMIR配准速度比SMIR快，但优势不明显，差距不大；在PET+PET1 配准时，RMIR配准速度比SMIR稍快，但优势不明显，差距不大；在PET+PET2 配准时，SMIR配准速度比RMIR稍快，但优势不明显，差距不大。因此，从 整体上来看，RMIR配准速度和SMIR相当，但图像较大时，RMIR稍逊SMIR。

根据表3-6和公式(9)和(10)，我们得到RMIR和SMIR误差，如表7所示。 从表7可知，RMIR和SMIR均能精确配准图像。RMIR(0.95)在配准CT+CT1和 PET+PET2时，精度最高；RMIR(1.1)在配准CT+CT2、MR+MR1和PET+PET1 时时，精度最高；SMIR只在配准MR+MR2时，精度最高。因此，RMIR达到 了SMIR的配准精度，甚至更高。因此，把RMIR作为相似性测度是一种有效 的方法。

表7单模态配准误差

下面进行多模态图像模拟试验，同样选取图2a、图2b和图2c作为参考 图像。实验分成三组：MR作为参考图像，CT1、CT2作为浮动图像，图像 大小均256×256，参考图22a、图22b，MR作为参考图像，PET1、PET2作 为浮动图像，图像大小均128×128；参考图23a、图23b，CT作为参考图像， PET1、PET2作为浮动图像，图像大小均128×128，在实验中也使用单纯形法作为多参数寻优方法，先使用RMIR进行多模态图像配准，然后使用SMIR 进行多模态图像配准。与单模态图像配准实验不同的是，在多模态图像配准 实验中，实际图像的配准空间变换参数的真实值是未知的，因此，通过计算 每幅实验图像的质心坐标和倾斜角，在此假定浮动图像和参考图像的质心坐 标和倾斜角之差的结果是相对准确的，并作为计算误差时的标准，如表8-9 所示。在以下配准结果图中，深灰色表示参考图像经过Canny算子提取边缘 的结果，浅灰色表示浮动图像经过Canny算子提取边缘的结果，白色部分表 示两幅图像配准后重叠部分。

表8实验图像质心坐标和倾斜角

表9多模态浮动图像变换参数

当α＝0.95时，将RMIR用于多模态配准得出实验结果4，实验结果如图 26a-f和表10所示。

表10实验4多模态配准性能

当α＝1.1时，将RMIR用于多模态配准得出实验结果5，实验结果如图 27a-f和表11所示。

表11实验5多模态配准性能

将SMIR用于多模态配准得出实验结果6，参考图28a-f和表12，

表12实验6基于互信息量多模态配准性能

从以上实验4、实验5和实验6结果来看，RMIR和SMIR相比较，在6个配 准实验中，在配准速度方面，SMIR、RMIR(0.95)和RMIR(1.1)三者非常接近， 差距均在0.5秒内。因此，从整体上来看，RMIR配准速度和SMIR相当。

根据表9-12和公式(9)和(10)，我们得到RMIR和SMIR的多模态配准方法 误差，如表13所示。

从表13可知，RMIR和SMIR均能准确配准图像。RMIR(0.95)在配准 MR+CT1和MR+CT2时，精度最高；RMIR(1.1)在配准PET+PET2时，精度最 高；SMIR在配准MR+PET1、MR+PET2和CT+PET1时，精度最高。因此， RMIR从整体上讲和SMIR有相似配准精度。因此，在多模态配准中，把RMI 作为相似性测度是一种同样有效的方法。

表13多模态配准误差

可以见得α的合适范围为：-0.1<α<0，0<α<1，1<α<1.3，可以获得 更好的配准性能，基于Renyi熵的医学图像匹配方法通过使用单纯形方法作 为寻优方法，可以减少算法运算量、提升配准速度，且具有较高的精度，并 且可以用于单模态图像配准，也可以用于多模态配准。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释 本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和 范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保 护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和 边界、或者这种范围和边界的等同形式内的全部变化和修改例。

Claims (4)

1.一种基于α-Renyi互信息的医学图像配准方法，其特征在于，包括如下步骤：

获取图像特征步骤，获取参考图像和浮动图像的图像特征；

获取图像灰度信息步骤，获取所述参考图像和所述浮动图像的灰度信息；

获取图像互信息量步骤，通过所述参考图像和所述浮动图像的所述图像特征和灰度信息获取Renyi互信息量，所述Renyi互信息量由下式获得：

RMI(R,F)＝H_α(R)+H_α(F)-H_α(R,F)

其中，RMI(R，F)为参考图像R和浮动图像F的Renyi互信息量，H_α(R)和H_α(F)分别为R和F的Renyi熵，H_α(R,F)是图像R和F的Renyi熵；

所述Renyi熵由下式获得：

其中，α为初始参数，p(k)为灰度值k出现的概率，L为灰度信息；

图像配准步骤，将所述获取图像互信息量步骤中获取的所述互信息量进行刚性配准，配准函数如下式：

其中，T(F)表示对浮动图像进行刚性变换；

通过所述配准函数获得最佳变换T_optimal，使得所述Renyi互信息达到最大值；

完成图像配准。

2.根据权利要求1所述的基于α-Renyi互信息的医学图像配准方法，其特征在于，

所述Renyi熵中参数α的取值范围为：-0.1<α<0。

3.根据权利要求1所述的基于α-Renyi互信息的医学图像配准方法，其特征在于，

所述Renyi熵中参数α的取值范围为：0<α<1。

4.根据权利要求1所述的基于α-Renyi互信息的医学图像配准方法，其特征在于，

所述Renyi熵中参数α的取值范围为：1<α<1.3。