CN102385748B

CN102385748B - 一种图像配准方法

Info

Publication number: CN102385748B
Application number: CN2010102701489A
Authority: CN
Inventors: 孙毅勇; 宫晶晶; 刘道志
Original assignee: Shanghai Microport Medical Group Co Ltd
Current assignee: Shanghai Microport EP MedTech Co Ltd
Priority date: 2010-08-31
Filing date: 2010-08-31
Publication date: 2013-12-25
Anticipated expiration: 2030-08-31
Also published as: ES2739052T3; WO2012028053A1; EP2613293A1; US9042614B2; US20130156281A1; CN102385748A; EP2613293B1; EP2613293A4

Abstract

本申请公开了一种图像配准方法，该方法包括：选取配准源图像{p_i}和配准目标图像{q_j}；对{p_i}施加随机扰动使其产生形变，得到{p_i’}，获取{q_j}上与{p_i’}上的点相对应的最近点的集合，即{q_j’}；根据预设的初始坐标变换H₀，对{p_i’}与{q_j’}进行迭代运算，得到{p_i’}与{q_j’}之间的坐标变换{H_l}，然后再计算与其相对应的平均距离{E_l}；判断平均距离{E_l}和预设的理想平均距离E_x的大小，且当{E_l}小于等于E_x时，配准终止。该方法无需手动调节参数，即可实现在配准源和配准目标在大范围初始偏移的情况下进行配准，而且能够有效地摆脱局部极限的束缚，提高了配准精度和配准成功率。

Description

一种图像配准方法

技术领域

本申请涉及图像处理技术领域，特别是涉及一种基于ICP算法的图像配准方法。

背景技术

随着新型传感器的不断涌现，人们获取图像的能力迅速提高，不同物理特性的传感器所产生的图像也不断增多。由于不同的图像传感器获取的图像数据存在明显的局限性和差异性，所以仅仅利用一种图像数据往往很难满足实际需求，为此需要通过图像融合技术将不同传感器获取的图像总和起来使用，以达到对图像中的目标更加全面清晰、准确的理解和认识。例如在医学上，为了通过对解剖和生理信息进行综合分析而通过对不同形式诸如计算机断层扫描成像（CT）、核磁共振成像（MRI）、超声（US）获得图像进行融合，以实现改进诊断过程。

图像配准技术是实现图像融合的重要前提，是图像融合首先要解决的问题。ICP（Iterative Closest Point，迭代最近点）算法是实现对图像进行配准的一种方法，ICP算法是一种基于自由形态曲面的配准方法，重复进行“确定对应点集，计算最优刚体变化”的过程，直到达到预设的配准收敛准则，最终坐标变化是每次变换的合成。但由于传统的ICP算法存在配准目标和配准源的初始偏移范围不能太大，可以收敛到某一局部极限，但通常不是全局最优（即无法挣脱局部极限）的缺点，因此图像配准的精度和成功率较低，无法满足实际需求。

通过对现有技术研究，申请人发现：现有的基于ICP算法的图像配准方法通常采用“给配准施加随机刚体扰动变换”和“给配准目标的三维坐标点施加随机扰动使其产生形变”两种方法达到挣脱传统ICP算法局部极限，使配准效果更好。但采用“给配准施加随机刚体扰动变换”这种方法需要在六维自由度空间采样，计算量大且费时；而采用“给配准目标的三维坐标点施加随机扰动使其产生形变”方法，其参数是在实际应用中根据不同的实验模型靠经验进行设定的，对于不同的实验模型，参数设定比较困难，设置错误就可能导致配准失败。

发明内容

有鉴于此，本申请实施例提供了一种图像配准方法，以实现无需手动调节参数，即可实现在配准源和配准目标在大范围初始偏移的情况下进行配准，而且能够有效地摆脱局部极限的束缚，提高了配准精度和配准成功率。

为实现上述目的，技术方案如下：

一种图像配准方法，包括：

选取配准源图像{p_i}和配准目标图像{q_j}，其中：0<i<M，0<j<N，M，N为正整数，且M<<N；

根据预设的随机扰动控制参数σ，对所述图像{p_i}施加随机扰动使其产生形变，得到图像{p_i’}，并获取所述图像{q_j}上与图像{p_i’}上的点相对应的最近点，得到的点的集合即为图像{q_j’}；

根据预设的初始坐标变换H₀，在预设步长L范围内，对图像{p_i’}与图像{q_j’}进行迭代运算，得到图像{p_i’}与图像{q_j’}之间的坐标变换{H_l}，0<l<L，根据坐标变换{H_l}计算与其相对应的平均距离{E_l}；

判断所述平均距离{E_l}和所述预设的理想平均距离E_x的大小，当所述平均距离{E_l}小于等于预设的理想平均距离E_x时，配准终止。

优选地，所述p_i和q_j为3*1的列向量。

优选地，所述获取所述图像{q_j}上与图像{p_i’}上的点相对应的最近点，具体为：

运用KD二叉树搜索所述图像{q_j}上与图像{p_i’}上的点相对应的最近点。

优选地，所述判断所述平均距离{E_l}和所述理想平均距离E_x的大小具体为：

依次判断所述平均距离{E_l}是否小于等于所述理想平均距离E_x。

优选地，所述判断所述平均距离{E_l}和所述理想平均距离E_x的大小具体体为：

选择所述平均距离{E_l}中的最小平均距离E_min；

判断所述最小平均距离E_min和所述理想平均距离E_x的大小。

优选地，当所述平均距离{E_l}大于预设的理想平均距离时，还包括：

依次判断所述平均距离{E_l}和初始平均距离E₀的大小，其中：所述初始平均距离E₀与所述初始坐标变换H₀相对应；

当所述平均距离{E_l}中任意一个值小于等于所述初始平均距离E₀，则终止所述平均距离{E_l}与所述初始平均距离E₀之间的判断，并将与当前平均距离E_g相对应的坐标变换H_g替换初始坐标变换H₀，且按照固定衰减速率将所述随机扰动控制参数σ进行衰减，根据坐标变换H_g和所述衰减后的随机扰动控制参数σ，重新计算平均距离{E_l’}，并判断所述平均距离{E_l’}是否小于等于所述预设的理想平均距离E_x的步骤；

当所述平均距离{E_l}中所有值均大于所述初始平均距离E₀时，则按照固定衰减速率将所述随机扰动控制参数σ衰减后，重新对所述图像{p_i}施加随机扰动，并利用初始坐标变换H₀，重新计算平均距离{E_l’’}，并判断所述平均距离{E_l’’}是否小于等于所述预设的理想平均距离E_x的步骤。

选择所述平均距离{E_l}中的最小平均距离E_min；

判断所述最小平均距离E_min和初始平均距离E₀的大小，其中：所述初始平均距离E₀与所述初始坐标变换H₀相对应；

当所述最小平均距离E_min小于等于所述初始平均距离E₀，则将与最小平均距离E_min相对应的坐标变换H_min替换初始坐标变换H₀，且按照固定衰减速率将所述随机扰动控制参数σ进行衰减，根据坐标变换H_min和所述衰减后的随机扰动控制参数σ，重新计算平均距离{E_l’}，并判断所述平均距离{E_l’}是否小于等于所述预设的理想平均距离E_x的步骤；

当所述最小平均距离E_min大于所述初始平均距离E₀时，则按照固定衰减速率将所述随机扰动控制参数σ衰减后，重新对所述图像{p_i}施加随机扰动，并利用初始坐标变换H₀，重新计算平均距离{E_l’’}，并判断所述平均距离{E_l’’}是否小于等于所述预设的理想平均距离E_x的步骤。

优选地，所述重新计算平均距离{E_l’}或重新计算平均距离{E_l’’}后，还包括：

判断迭代循环次数是否大于预设的迭代循环次数K，并且当所述迭代循环次数大于预设的迭代循环次数K时，则配准终止，并报错。

优选地，所述随机扰动控制参数σ的固定衰减速率为

。

优选地，所述步长L的范围为50～200。

由以上本申请实施例提供的技术方案可见，该方法中的随机扰动控制参数σ在整个迭代过程中可以自动衰减。当预设的初始变换是最佳起始位置时，无需担心因为增加了随机扰动而将模型变换到更差的位置；当初始变换本身不是最佳时，还可以通过增加随机扰动使模型挣脱局部极值的束缚而变换到更好的位置进行迭代，即通过搜索可以接近全局最优。同时，该方法中随机扰动控制参数σ不是根据经验设定的，而是考虑了全局搜索性能与迭代循环速度后综合设定的。对于多组心腔模型，随机扰动控制参数σ可以设置为同一值，且都得到较好的配准效果。

因此，本申请实施例提供的图像配准方法，在无需手动调节参数的情况下，不仅可以实现配准源和配准目标在大范围初始偏移的情况下进行配准，而且能够有效地摆脱局部极限的束缚，提高了配准精度和配准成功率。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本申请实施例一提供的一种图像配准方法的流程示意图；

图2为本申请实施例二提供的一种图像配准方法的流程示意图；

图3为本申请实施例三提供的一种图像配准方法的流程示意图；

图4为本申请实施例四提供的一种图像配准方法的流程示意图；

图5为本申请实施例四提供的另一种图像配准方法的流程示意图。

具体实施方式

为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。

实施例一：

图1为本申请实施例一提供的一种图像配准方法的流程示意图。

如图1所示，该方法包括：

步骤S101：选取配准源图像{p_i}和配准目标图像{q_j}。

配准源图像为{p_i|i=1,……,M}，配准目标图像为{q_j|j=1,……,N}，M<<N，配准源图像{p_i}和配准目标图像{q_j}处于不同的坐标系，

且{p_i}=H*{q_j}，

H = [\begin{matrix} R & 0 \\ T & 1 \end{matrix}],

其中：0<i<M，0<j<N，M，N为正整数，且M<<N；R为3*3的旋转矩阵，T为3*1的平移向量；

步骤S102：根据预设的随机扰动控制参数σ，对图像{p_i}施加随机扰动使其产生形变，得到图像{p_i’}，并获取图像{q_j}上与图像{p_i’}上的点相对应的最近点，得到的点的集合即为图像{q_j’}。

首先对配准源图像{p_i}施加随机扰动，具体为：

给配准源图像{p_i}中每个点增加满足（0，σ²）正态分布的随机扰动n_i，使得配准源图像{p_i}产生形变，

即p_i’=p_i+n_i，

其中，n_i为3*1的列向量。

然后，运用KD二叉树搜索配准目标图像{q_j|j=1,……,N}与形变后的配准源图像{p_i’|i=1,……,M}上相对应的最近点，搜索得到的点的集合即为{q_j’|j=1,……,N}。

并根据

计算3*3矩阵X，并根据X=UDV^t对矩阵X进行奇异值分解，其中：

{p_{i}}^{''} = {p_{i}}^{'} - \overset{&OverBar;}{p},

{q_{i}}^{''} = {q_{j}}^{'} - \overset{&OverBar;}{q},

\overset{&OverBar;}{p} = \frac{1}{M} Σ_{i = 1}^{M} {p_{i}}^{'},

\overset{&OverBar;}{q} = \frac{1}{M} Σ_{j = 1}^{M} {q_{j}}^{'},

U和V是3*3的正交矩阵，D为由奇异值组成的对角矩阵。

最后，分别计算旋转矩阵R和平移向量T，其中：R=VU^t，

步骤S103：根据预设的初始坐标变换H₀，在预设步长L范围内，对图像{p_i’}与图像{q_j’}进行迭代运算，得到图像{p_i’}与图像{q_j’}之间的坐标变换{H_l}，0<l<L，根据坐标变换{H_l}计算与其相对应的平均距离{E_l}。

根据{p_i}=H*{q_j}，且

H = [\begin{matrix} R & 0 \\ T & 1 \end{matrix}],

计算{p_i}和{q_j}间的坐标变换H，并根据坐标变换H计算{p_i}和{q_j}间的平均距离E，其中：

E = \frac{1}{M} Σ_{i, j = 1}^{M} {| | {p_{i}}^{'} - ({Rq}_{j}^{'} + T) | |}^{2} .

并且对图像{p_i’}与图像{q_j’}进行L次迭代运算，并记录L次迭代过程中图像{p_i}和图像{q_j}之间的坐标变换{H_l|l=1,……,L}，并根据坐标变换{H_l|l=1,……,L}计算得到图像{p_i}和图像{q_j}之间平均距离{E_l|l=1,……,L}，其中：0<l<L。

步骤104：判断平均距离{E_l}是否大于预设的理想平均距离E_x，如果否，配准终止。

如果平均距离{E_l}的中任意一个值小于等于预设的理想平均距离E_x，则平均距离{E_l}的中的这个值即为理想坐标变换，所以配准可以终止。

实施例二：

图2为本申请实施例二提供的一种图像配准方法的流程示意图。

如图2所示，本申请实施例中步骤S201～步骤S204与实施例一中的步骤S101～步骤S104相同，再次不再赘述，在步骤S204中如果平均距离{E_l}大于预设的理想平均距离E_x，则该方法还可以包括：

步骤S205：依次判断平均距离{E_l}是否大于初始平均距离E₀，其中：初始平均距离E₀与初始坐标变换H₀相对应；

步骤S206：如果否，则终止平均距离{E_l}与初始平均距离E₀之间的判断，并将与当前平均距离E_g相对应的坐标变换H_g替换初始坐标变换H₀，且按照固定衰减速率将随机扰动控制参数σ进行衰减，根据坐标变换H_g和衰减后的随机扰动控制参数σ，重新计算平均距离{E_l’}，并判断平均距离{E_l’}是否小于等于预设的理想平均距离E_x的步骤。

步骤S207：如果是，则按照固定衰减速率将随机扰动控制参数σ衰减后，重新对图像{p_i}施加随机扰动，并利用初始坐标变换H₀，重新计算平均距离{E_l’’}，并判断平均距离{E_l’’}是否小于等于预设的理想平均距离E_x的步骤。

算法的全局搜索性能与迭代速度呈反关系通过设置随机扰动控制参数σ的范围和设置步长L可以平衡全局搜索性能和迭代速度之间的关系，在本申请实施例中，预设随机扰动控制参数σ的为配准源图像{p_i}或者配准目标图像{q_j}几何比例的10%。步长L的范围为50～200。

实施例三：

图3为本申请实施例三提供的一种图像配准方法的流程示意图。

如图3所示，本申请实施例中步骤S301～步骤S304与实施例一中的步骤S101～步骤S104相同，再次不再赘述，在步骤S304中判断平均距离{E_l}和预设的理想平均距离E_x的大小之后，该方法还可以包括：

步骤S305：选择平均距离{E_l}中的最小平均距离E_min；

步骤S306：判断最小平均距离E_min是否大于初始平均距离E₀，其中：初始平均距离E₀与初始坐标变换H₀相对应；

步骤S307：如果否，则将与最小平均距离E_min相对应的坐标变换H_min替换初始坐标变换H₀，且按照固定衰减速率将随机扰动控制参数σ进行衰减，根据坐标变换H_min和衰减后的随机扰动控制参数σ，重新计算平均距离{E_l’}，并判断平均距离{E_l’}是否小于等于预设的理想平均距离E_x的步骤；

步骤S308：如果是，则按照固定衰减速率将随机扰动控制参数σ衰减后，重新对图像{p_i}施加随机扰动，并利用初始坐标变换H₀，重新计算平均距离{E_l’’}，并判断平均距离{E_l’’}是否小于等于预设的理想平均距离E_x的步骤。

步骤S305和步骤S205所不同的只是如何判断平均距离{E_l}与初始平均距离E₀的大小，步骤S205中，将平均距离{E_l}依次与初始平均距离E₀进行比较，而步骤S305中，则是首先选取平均距离{E_l}中的最小平均距离E_min，然后用最小平均距离E_min与初始平均距离E₀进行比较。

实施例四：

图4为本申请实施例四提供的一种图像配准方法的流程示意图。图5为本申请实施例四提供的另一种图像配准方法的流程示意图。如图4和图5所示，在实施例二中步骤S206或S207之后，或者在实施例三步骤S307或S308之后，该方法还可以包括：

步骤S400：判断迭代循环次数是否大于预设的迭代循环次数K，如果是，则配准终止，并报错。

预设迭代循环次数K，是为了将迭代循环控制在合理的范围内，避免无法配准的情况下，浪费大量的时间进行运算。在本申请实施例中，预设的迭代循环次数K选择100。

由以上本申请实施例提供的技术方案可见，该方法中的随机扰动控制参数σ在整个迭代过程中可以自动衰减，当预设的初始变换是最佳起始位置时，无需担心因为增加了随机扰动而将模型变换到更差的位置；并且当初始变换本身不是最佳时，还可以通过增加随机扰动可以使模型挣脱局部极值的束缚而变换到更好的位置进行迭代，即可以搜索可以接近全局最优。同时，该方法中随机扰动控制参数σ不是根据经验设定的，而是考虑了全局搜索性能与迭代循环速度后综合设定的，对于多组心腔模型，随机扰动控制参数σ可以设置为同一值，且都得到较好的配准效果。

因此，本申请实施例提供的图像配准方法，在无需手动调节参数的情况下，不仅可以实现在配准源和配准目标在大范围初始偏移的情况下进行配准，而且能够有效地摆脱局部极限的束缚，提高了配准精度和配准成功率。

以上仅是本申请的优选实施方式，使本领域技术人员能够理解或实现本申请。对这些实施例的多种修改对本领域的技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本申请的精神或范围的情况下，在其它实施例中实现。因此，本申请将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims

1.一种图像配准方法，其特征在于，包括：

选取配准源图像{p_i}和配准目标图像{q_j}，且{p_i}=H*{q_j}，

H = [\begin{matrix} R & 0 \\ T & 1 \end{matrix}],

其中：0<i<M，0<j<N，M，N为正整数，且M<<N，R为3*3的旋转矩阵，T为3*1的平移向量；

根据预设的初始坐标变换H₀，在预设步长L范围内，对图像{p_i’}与图像{q_j’}进行迭代运算，得到图像{p_i’}与图像{q_j’}之间的坐标变换{H_l}，0<l<L，根据坐标变换{H_l}计算与其相对应的平均距离{E_l}，其中，所述平均距离为：

E_{1} = \frac{1}{M} Σ_{i, j = 1}^{M} {| | {p_{i}}^{'} - ({Rq}_{j}^{'} + T) | |}^{2};

2.根据权利要求1所述的方法，其特征在于，所述p_i和q_j为3*1的列向量。

3.根据权利要求2所述的方法，其特征在于，所述获取所述图像{q_j}上与图像{p_i’}上的点相对应的最近点，具体为：

4.根据权利要求3所述的方法，其特征在于，所述判断所述平均距离{E_l}和所述理想平均距离E_x的大小具体为：

5.根据权利要求3所述的方法，其特征在于，所述判断所述平均距离{E_l}和所述理想平均距离E_x的大小具体为：

选择所述平均距离{E_l}中的最小平均距离E_min；

判断所述最小平均距离E_min和所述理想平均距离E_x的大小。

6.根据权利要求1所述的方法，其特征在于，当所述平均距离{E_l}大于预设的理想平均距离时，还包括：

7.根据权利要求1所述的方法，其特征在于，当所述平均距离{E_l}大于预设的理想平均距离时，还包括：

选择所述平均距离{E_l}中的最小平均距离E_min；

8.根据权利要求6或7所述的方法，其特征在于，所述重新计算平均距离{E_l’}或重新计算平均距离{E_l’’}后，还包括：

9.根据权利要求6或7所述的方法，其特征在于，所述随机扰动控制参数σ的固定衰减速率为

。

10.根据权利要求1-7任一项所述的方法，其特征在于，所述步长L的范围为50～200。