CN1920882A - 基于显著区域特征的医学图像三维多模配准的***和方法 - Google Patents

Abstract

一种用于对准一对图像的方法包括：提供一对图像；在第一图像和第二图像中标识(60)显著特征区域，其中每个区域与空间尺度相关；通过每个区域的中心点来表示(64)特征区域；根据局部强度将一幅图像的特征点与另一幅图像的特征点进行配准(65)；通过相似性度量对所述特征对进行排序；并且通过将中心点细化(68)到子像素精度来优化特征对的联合对应集合。

Description

基于显著区域特征的医学图像 三维多模配准的***和方法

技术领域

本发明涉及三维数字医学图像的配准技术。

背景技术

本申请要求了Xu等人于2005年8月24日提交的发明名称为“Method and System for Robust Salient Region BasedRegistration of 3D Medical Images(基于稳定的显著区域的三维医学图像配准的方法和***)”的美国临时申请No.60/710,834的优先权，该美国临时申请的内容在此被引入作为参考。

在医学图像处理中，配准已经成为在两幅或者多幅图像的空间定位之间产生映射的基本任务，并且能够被用在各种应用中。对准变换的主要要求是最优覆盖相应的图像内容。现有的配准方法能够被分类为基于特征(例如，标志点(landmark)、边缘、标记)的配准方法、基于强度的配准方法、或结合上述两种方法的多个方面的混合方法。基于标志点的方法能够导致在远离特征的区域中的粗劣的对准。仅仅基于图像强度的方法由于优化过程的特性易陷入到决非理想解决方案的局部最优中。混合技术使用多种这样的属性的组合并且尤其优选于来自不同模态的图像配准，其中图像强度或几何信息单独并没有提供准确的测量基础。例如，和附加信息通道配对的互信息可以改善MR和PET图像的配准结果，所述附加信息通道由区域标记信息组成。混合配准技术例如在低对比度的腹部区域、血浆凝胶电泳或蛋白质成像等中是公知的。其它应用领域例如是创建图谱或标准数据库，这些图谱或标准数据库适于图像或对象分析、允许医生获取对疾病进展的了解的同一患者或患者间研究、或在癌症治疗期间基于时间的追踪研究。

对同一个受检者使用不同的成像***能够获取更多的信息，但是另一方面这要求用于合适解释的多模配准技术。补充信息的增加被各种医学成像***所利用，这些医学成像***大致可被分成两大类：提取形态学信息的解剖学成像(例如，X射线、计算机断层扫描(CT)、磁共振成像(MRI)、超声波(US))，和可视化下层解剖结构的新陈代谢信息的功能成像(例如，单光子发射计算机断层扫描(SPECT)、正电子发射断层扫描(PET)、功能MRI(fMRI))。在多模图像配准中，不同类型图像的组合对医生是有利的。例如，CT图像的特征在于良好的空间分辨率，而PET图像描述了下层组织的功能。因此，通过与部分缺乏空间分辨率的相应的PET图像融合，CT图像中功能信息的不足能够被弥补。

发明内容

如在此所描述的本发明的示例性实施方式通常包括用于配准方法的方法和***，上述配准方法从两幅图像中自动提取三维区域，在这两幅图像之间找到相应对，并建立刚性配准变换，以便使医学应用中的融合结果可视化。那些固有特征的尺度、平移和旋转不变性适于三维，以估计下层的单模或多模三维医学图像之间的变换。根据本发明实施方式的方法结合了基于特征和基于强度方法的有利方面，并且包括每一幅图像上的三维显著区域特征集合的自动提取、相对应集合的估计及其子像素的准确细化，上述细化包括排除逸出值(outlier)。基于区域生长的方法被用于提取三维显著区域特征，这就减少了特征聚类和相对应的搜索空间的复杂度。

根据本发明实施方式的算法以通过使用kD树结构进行显著区域的快速聚类为特征，并且该算法使用局部强度驱动的三维显著特征区域的配准来提高优化。根据本发明实施方式的算法的附加特征包括在两对显著性区域的质心上使用迭代最近点(iterative closest point)算法来执行初始姿态估计和使用基于局部强度的配准进行质心的局部细化以获得子像素精度。

根据本发明实施方式的算法是全自动的、稳定的、多模刚性图像配准算法，该算法能够在任意姿态下成功配准三维图像。在几何构型约束条件的约束下通过各种尺度不变的三维特征对图像建模。完全利用同一图像上的特征之间的相对构型约束条件来继续进行多个三维显著特征对之间的联合对应。与单独的特征对之间的一致性相比，通过联合对应施加的严格的几何约束使得极不可能发生误匹配。当联合对应逐渐增加直到全局图像对准质量不能通过增加新的对进一步提高时，从联合对应集合中估计得到的变换收敛于全局最优。这通过合适的收敛标准来获得。

根据本发明的方面，提供了一种用于对准一对图像的方法，该方法包括如下步骤：提供具有第一图像和第二图像的一对图像，其中所述图像包括多个对应于三维空间中的像素域的强度；在第一图像和第二图像中均标识显著特征区域，其中每个区域都与空间尺度相关；通过每个区域的中心点代表特征区域；根据局部强度将一幅图像的特征点与另一幅图像的特征点进行配准；通过相似性度量来对所述特征对进行排序；通过将中心点细化到子像素精度来优化特征对的联合对应集合。

根据本发明的另一方面，该方法包括：在kD树中表示一幅图像的显著特征区域中心点；查询所述kD树的每个特征来找到一组最邻近的特征；以及从所述树中移除具有较低显著性值的那些最邻近特征和在所述每个特征的尺度内具有中心点的那些最邻近特征，其中获得显著特征区域在所述图像中的基本上均匀的分布。

根据本发明的另一方面，所述空间尺度是包括所述特征区域的球形的半径。

根据本发明的另一方面，所述kD树使用所述显著特征区域中心点的图像像素索引作为树叶，并且其中从特征区域到最邻近特征区域的距离以图像索引为单位。

根据本发明的另一方面，根据局部强度对特征点进行配准还包括：使用所述第一图像和所述第二图像之间的迭代最近点变换来估计初始配准；将所述第二图像的所有特征变换到所述第一图像的坐标空间；将所述所变换的特征存储在kD树中，并且根据预定的选择标准查询所述树的所述第一图像中的每个特征以选择所述第二图像中的那些最邻近特征；并且在所述第一图像和所述第二图像中测试所述特征的所选择的特征对的平移不变性、旋转不变性和全局图像相似性度量，其中通过所述所选择的特征对的全局图像相似性度量值来对所选择的特征对进行排序。

根据本发明的另一方面，所述迭代最近点变换最小化每组特征点之间的均方误差；

根据本发明的另一方面，测试所述平移不变性包括估计 ${\widehat{\mathrm{\Θ}}}_{i,j}^T=p_i-p_j,$ 其中p_i和p_j是第i个第一图像特征和第j个第二图像特征在物理空间的中心位置坐标。

根据本发明的另一方面，测试所述旋转不变性包括估计

${\widehat{\mathrm{\Θ}}}_{i,j}^R=\underset{{\mathrm{\Θ}}^R}{\arg \max }\mathrm{ECC}(f_i,f_j^{T_{{\mathrm{\Θ}}^R}}),$

其中(f_i，f_j)分别表示所述第一图像和所述第二图像中的所述特征对，

$\mathrm{ECC}(f_i,f_j)=2-\frac{2H(f_i,f_j)}{H\left(f_i\right)+H\left(f_j\right)},$

其中H表示相对于在体素位置x和所述空间尺度s周围的球形邻近特征区域f_s内的图像强度值的熵，H被定义为

$H_D(s,x)=-{\∫}_{R_s}p(i,s,x){\log }_{2}p(i,s,x)\mathrm{di},$

其中p(i，s，x)是被包含在f中的图像强度值i的概率密度函数，其中H(f_i，f_j)是联合微分熵，H(f_i，f_j)被定义为：

$H(f_i,f_j)=-{\∫}_{f_i,f_j}p(f_i,f_j){\log }_{2}p(f_i,f_j)\mathrm{dIdJ},$

其中p(f_i，f_j)是特征区域f_i和f_j中的图像强度的联合概率密度，并且I和J分别在所述第一和第二图像中的可能的强度值的集合中取值。

根据本发明的另一方面，测试所述全局图像相似性度量包括估计

$L_{\mathrm{global}}\left(c_{i,j}\right)=\mathrm{ECC}(I_r,I_t^{T_{{\widehat{\mathrm{\Θ}}}_{i,j}}}),$

其中I_r表示所述第一图像，

表示所述第二图像到所述第一图像的坐标空间上的所述变换。

$\mathrm{ECC}(I_r,I_t^{T_{{\widehat{\mathrm{\Θ}}}_{i,j}}})=2-\frac{2H(I_r,I_t^{T_{{\widehat{\mathrm{\Θ}}}_{i,j}}})}{H\left(I_r\right)+H\left(I_t^{T_{{\widehat{\mathrm{\Θ}}}_{i,j}}}\right)},$

其中H表示相对于在体素位置x和所述空间尺度s周围的所述图像中的一幅图像内的图像强度值的熵，H被定义为：

H(s，x)＝-∫_Ip(i，s，x)log₂p(i，s，x)di，

其中p(i，s，x)是被包含在I中的图像强度值i的概率密度函数，其中是联合微分熵，

被定义为：

$H(I_r,I_t^{T_{{\widehat{\mathrm{\Θ}}}_{i,j}}})=-{\∫}_{I_r,I_t^{T_{{\widehat{\mathrm{\Θ}}}_{i,j}}}}p(I_r,I_t^{T_{{\widehat{\mathrm{\Θ}}}_{i,j}}}){\log }_{2}p(I_r,I_t^{T_{{\widehat{\mathrm{\Θ}}}_{i,j}}})\mathrm{dIdJ},$

其中 是图像I_r和 中的图像强度的联合概率密度，并且I和J分别在所述第一和第二图像中的可能的强度值的集合中取值，并且其中L_global在所述第一和第二图像的整个重叠域上被评估。

根据本发明的另一方面，优化特征对的联合对应集合还包括：利用根据所述相似性度量为最相似的特征对来对所述联合对应集合进行初始化；针对所述联合对应集合与没有已被包括在所述联合对应集合中的每个特征对的并集(union)来估计所述相似性度量；选择最大化所述并集的相似性度量的特征对，其中如果所述最大化特征对与所述联合对应集合的并集的相似性度量大于该联合对应集合的相似性度量，则所述最大化特征对利用局部刚性变换以子像素精度进行配准并且被增加到所述联合对应集合中。

根据本发明的另一方面，通过使用迭代最近点过程计算特征对之间的配准变换来使相似性度量最大化。

根据本发明的另一方面，如果所述最大化特征对与所述联合对应集合的并集的相似性度量小于或等于该联合对应集合的相似性度量，则提供配准变换，所述配准变换由最大化所述相似性度量的特征对之间的配准变换计算得到。

根据本发明的另一方面，提供一种计算机可读的程序存储装置，该计算机可读的程序存储装置确实地包含了可由计算机执行的程序指令，以执行用于对准一对图像的方法步骤。

附图说明

图1(a)-(e)根据本发明的实施方式示出显著特征区域的示例性集合。

图2(a)-(d)根据本发明的实施方式示出CT体积的切片(左部)，上述切片已经被平移、旋转并且被覆盖到原始切片上。

图3根据本发明的实施方式示出用于优化联合对应集合和配准变换的EM型算法的流程图。

图4(a)-(c)根据本发明的实施方式是显示所有测量距离和标准偏差的表格，所述测量距离分别针对PET-CT、CT-CT和SPECT-CT体积对在x、y和z方向上以cm为单位给出。

图5(a)-(c)示出使用根据本发明实施方式的算法获得的来自融合配准结果图像的三个切片，上述融合配准结果图像分别来自PET-CT图像对、具有强度伪像的CT-CT图像对和SPECT-CT图像对。

图6是根据本发明实施方式的三维配准过程的流程图。

图7是用于实现根据本发明实施方式的三维配准过程的示例性计算机***的框图。

具体实施方式

如在此所描述的本发明的示例性实施方式通常包括用于自动提取的尺度不变、旋转不变和平移不变的三维显著区域特征的构型匹配的***和方法。因此，虽然本发明容许各种修改和替换形式，但是仍通过附图中的实例示出其具体实施方式并将在此详细描述该具体实施方式。但是，应当理解，没有意图将本发明限于所公开的特定形式，而是相反，本发明要涵盖落入本发明的精神和范围内的所有修改、等价物和替换方式。

如在此所使用的那样，术语“图像”是指由离散图像元素(例如，二维图像的像素和三维图像的体素)组成的多维数据。例如，图像可以是由计算机断层扫描、磁共振成像、超声波或本领域技术人员所熟知的任何其它医学成像***所采集的受检者的医学图像。图像也可以从非医学环境中提供，诸如由遥感***、电子显微技术来提供等。尽管图像能够被看作是从R³到R的函数，但是本发明方法没有被限制到这样的图像上，而是能够被应用到任意维数的图像(例如，二维图片或三维体积)。对于二维或三维图像，图像域通常是二维或三维的矩形阵列，其中每个像素或体素能够关于2个或3个相互正交轴的集合来编址。如在此所使用的术语“数字”和“数字化”指(如适当)通过数字采集***或通过模拟图像转换所获取的数字或数字化格式的图像或体积。

根据本发明的实施方式，显著性描述(saliency description)被用于从三维图像中自动提取特征。在图6中呈现概述了根据本发明实施方式的三维配准算法的流程图。下文会对其中的步骤的细节进行详细描述。现在参考该图，在步骤60，通过形成球形区域，三维显著区域特征被检测到。所得到的每个显著区域特征均具有下述属性：(1)区域中心，(2)区域尺度(半径)和(3)区域的显著性值。

在步骤61，通过将一个集合的特征区域中心点存储到允许快速查询空间位置的kD树结构中，获得显著区域特征在图像空间上的均匀分布。在步骤62，针对每个特征来对树进行查询，以找到最邻近特征。因为是通过其尺度(半径)对区域进行描述，所以在步骤63中移除具有较低显著性的那些特征和位于当前特征的尺度内的中心点。所得到的聚类显著特征的子集被从整个集合中移除。没有显著区域特征位于具有较高显著性值的显著区域特征的尺度内。这种方法避免了能够在图像区域中出现的显著区域特征的聚类，上述图像区域具有大且是局部最大的显著性值。因此，根据那些特征的随后的图像配准可能出现偏差并导致错误解。

旋转不变性的参数通过三维显著特征区域的局部的、强度驱动的配准进行估计。局部变换的参数搜索空间被限制到旋转参数上。在这些通过测试每个组合来发现显著区域特征之间的一致性具有高计算复杂度，这些显著区域特征被包含在两个不相关的集合中。在两个三维显著区域特征集合的情况下，在步骤64，这些特征被缩减到其中心点。在步骤65，利用迭代最近点(ICP)算法对所得到的两个点云进行配准，该迭代最近点(ICP)算法最小化集合之间的均方误差。在步骤66，所得到的变换被用于将对准的图像特征集合映射到参考图像空间中。在步骤67，所变换的对准特征能够被存储到kD树结构中，以便允许快速查询相对于参考图像显著区域特征的最近N个被变换的对准特征。根据初始ICP变换产生点集之间的稳定映射的假设，对所假定的一致性的测试能够替代对准图像显著特征区域的整个集合而被缩减到每个参考特征的N个相邻特征。

为了在显著特征区域一致性之间实现更大的精度，在步骤68，通过显著特征区域之间的局部强度驱动的配准来以子像素精度对准这些特征中心。在这种情况下，针对预期最大化优化过程的参数搜索空间包括平移、旋转和尺度参数。

显著区域特征的原理是大量局部不可预测性或信号复杂度相对于某个尺度的表达，其中尺度是指体素周围的球形区域的半径。这种方法根据不同尺度的圆形区域的香农(Shannon)熵来区分感兴趣点，并且假设来自不同的相应解剖学结构或功能结构的体素具有相似的显著性值。显著性描述的局部属性为图像配准提供了下述主要优点：即使图像并没有重叠，在不同图像之间相对应的显著区域特征相对于总空间变换是不变的。

与相符的CT图像相比，SPECT和PET图像的观察表示，局部显著性最大值的位置经常可以在相应的感兴趣结构内被局部平移。这能够通过局部刚性配准步骤来解决，该局部刚性配准步骤包括在图像内容没有发生变形的情况下根据局部强度的相关性将区域中心子像素准确地平移到相对应的位置。这个步骤保持了针对上面所提及的相对应特征的类似显著值的基本假设。

针对图像强度范围D，如下定义显著性：

A_D(s_p，x)＝H_D(s_p，x)·W_D(s_p，x)，

其中H_D表示相对于在体素位置x和尺度(半径)s周围的球形邻近区域R_s内的图像强度值i∈D的熵：

$H_D(s,x)=-{\∫}_{R_s}p(i,s,x){\log }_{2}p(i,s,x)\mathrm{di}.---\left(1\right)$

这里，p(i，s，x)是针对被包含在R_s中的图像强度值的描述符i的概率密度函数(PDF)。W_D(s，x)是PDF之间的相似性相对于尺度的度量。W_D(s，x)随着PDF的不相似性的增大而增大：

$W_D(s,x)=s{\∫}_{R_s}\left|\frac{\∂}{\∂s}p(i,s,x)\right|\mathrm{di}.$

在x处形成H_D的局部峰值的尺度s_p由下式给出：

$s_p=\{s:\frac{\&PartialD;H_D(s,x)}{\&PartialD;s}=0,\frac{{\&PartialD;}^2{H}_D(s,x)}{{\&PartialD;}^{2}s}<0\}.---\left(2\right)$

在针对每个体素求解方程式(1)之后，结果是与输入图像大小相等的两个临时图像，这两个临时图像必须被分析：一个图像包含实际的显著性值，而另一个图像包括来自方程式(2)的尺度值。利用区域生长搜索方法，可以从显著性图像中提取局部最优和最具描述性的显著区域点。对于搜索空间的减小，全局显著性阈值δ被用作下限。根据本发明的实施方式，将全局显著性阈值δ设成平均显著性的一半 $(\mathrm{\&delta;}=\frac{1}{2}\stackrel{\&OverBar;}{A})$ 的经验设置对于排除无意义区域取得了良好结果。

基于区域特征的kD树结构化的最近点算法能够确定局部显著性最大值的位置，这形成了体素位置的列表，这些体素位置根据其显著性值进行排序。这种方法避免了局部最大值的聚类，如果全局阈值被应用并且仅根据降序的显著性值提取特征，则例如出现局部最大值的聚类。使用所提取的局部显著性最大值的区域中心的索引作为树叶，创建kD树。通过查询特征的区域中心索引，特定特征的K个最近邻能够被有效地找到。因此，到所返回的特征的距离不是以物理单位为量纲，而是以图像索引为量纲。尺度参数能够被用作最小距离要求：距离等于或小于所查询的特征的尺度和具有较低显著性的所有所返回的特征被从特征集中移除。这种限制能够被应用到整个集合中，以便移除聚类区域。如果结果集合要求特定的大小，那么可以在列表中添加满足距离标准的、有较低显著性的特征。如果特征中心没有位于具有较高显著性值的特征的区域内，那么该特征被保留在该集合中。所得到的三维显著区域特征的集合均匀分布，并且为随后的特征对应搜索提供了情况良好的初始集合。

图1(a)-(e)示出显著特征区域的示例性集合。显著特征区域在图1(a)-(b)中用白色圆圈表示。图1(a)示出对显著区域特征的集合进行聚类的效果，同时图1(b)示出通过根据本发明的实施方式的方法所选择的显著区域特征。图1(c)、1(d)和1(e)示出如在分别从CT图像、PET图像和MR图像中提取之后被可视化的最重要的三维显著区域特征。显著区域特征在图1(c)-(e)中表现为球形泡。利用特定传递函数来限幅该体积，以便在三维中对特征的位置进行可视化，而提取自身已经在整个强度范围上被执行。

三维配准的下一步(被称作区域分量匹配步骤)估计所假设的两幅图像的特征之间的对应的集合。假设I_r是指参考图像且I_t是指模板图像，假设N_r是从I_r中所提取的特征的数目且N_t是从I_t中所提取的特征的数目。所有假设特征对应的集合是C＝{c_i，j}，其中i∈[1，N_r]，j∈[1，N_t]，|C|＝N_r×N_t并且其中c_i，j＝(f_i，f_j)是一对I_r中的特征f_i和I_t中的特征f_j。

参数集合Θ限定了对准两幅图像的变换T，并且能够根据f_i和f_j之间的平移、尺度和旋转不变性进行估计。f_i和f_j之间的平移部分能够根据下式直接估计： ${\widehat{\mathrm{\&Theta;}}}_{i,j}^T=p_i-p_j,$ 其中p_i和p_j是物理空间中的第i个参考特征和第j个模板特征的中心位置。这种情况下的尺度不变性不是必需的，因为对于根据本发明实施方式的三维医学图像，体素大小被提供在DICOM(医学数字成像和通信)的报头中。为了实现旋转不变性，根据三维显著特征区域的强度值，通过这些三维显著特征区域的局部刚体配准来估计旋转参数。优化被限制在旋转参数子空间Θ^R中，并且由强度相似性度量来驱动，熵校正系数作为标准化互信息的特定形式由下式确定：

$\mathrm{ECC}(A,B)=2-\frac{2H(A,B)}{H\left(A\right)+H\left(B\right)},$

其中联合微分熵H(A，B)能够被定义成：

H(A，B)＝-∫_A，Bp(A，B)log₂p(A，B)dIdJ，

其中积分域是在区域R_i和R_j ^Θ上，p(A，B)是图像强度在区域A和B中的联合概率密度，而I和J分别采用I_f和I_m中的可能强度值集合中的值。这个系数为重叠域提供了提高的稳定性，以及一些另外的有利属性：增大的值表示图像之间增大的相关性，反之亦然即减小的值表示图像之间减小的相关性。因此，旋转不变性能够用公式表示为Θ^R的优化问题：

${\widehat{\mathrm{\&Theta;}}}_{i,j}^R=\underset{{\mathrm{\&Theta;}}^R}{\arg \max }\mathrm{ECC}(f_i,f_j^{T_{{\mathrm{\&Theta;}}^R}}).$

全局图像相似性度量L_global被用来估计M对中的每对的质量：

$L_{\mathrm{global}}\left(c_{i,j}\right)=\mathrm{ECC}(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}),$

其中L_global在两个图像的整个重叠域上被估计，而不是正好在局部特征区域上被估计。

隔开大空间距离的特征对不可能是相对应的并且能够从所假设的相对应的集合中被移除，这就减少了联合对应的搜索空间。因此，对应搜索空间能够从所有对的组合中被减少到仅在局部最邻近特征之间的组合。通过将集合看作是区域中心位置周围的点云来计算参考区域特征集合与模板区域特征集合之间的ICP变换，相邻集合能够被估计。这种ICP算法利用局部最小均方误差(MSE)来对准特征集合。该结果被用来将所有模板特征变换成参考图像的坐标空间，并将所有模板特征存储到新的kD树中。然后，针对参考图像中的每个显著特征，大致最邻近的特征能够在树的快速搜索中被确定。所变换的模板特征的最邻近特征的数目N_n是比集合的整个基数小得多的值：N_n＜＜N_t，所述最邻近特征和每个参考特征组合在一起。根据初始ICP变换是实际对准变换的良好近似的假设，并且根据如果特征相隔更大距离那么这些特征更不可能相对应的假设，这将复杂度减小到N_r×N_n。N_r×N_n个特征对在其平移不变性

旋转不变性 和全局图像相似性度量L_global(c_i，j)的基础上针对其对应质量进行测试，所述平移不变性 旋转不变性和全局图像相似性度量L_global(c_i，j)能够通过在整个对准的图像上施加局部特征变换来实现。

在根据本发明的实施方式执行的实验中，N_n＝(1/10)％N_t的邻域大小已经被成功地用来建立联合对应的初始搜索空间。此外，假定对应通过全局相似性度量L_global来排序，这在所估计的对应集合中导致更少的逸出值。

图2(a)-(d)示出CT体积的切片(左部)，该切片已经被平移、旋转和覆盖到原始切片上。圆圈表示具有特定尺度的显著特征区域。根据L_global已经执行了排序，并且整个集合包含很少的逸出值。为了清楚起见，针对每种方法只示出前四个对应。

大小为M的假定对应 $C=\{c_{i_1,j_1},...,c_{i_M,j_M}\}$ 的集合被用来估计两幅图像之间的变换T，所述集合以算法的特征对应搜索计算出来。这种变换不是十分准确，因为其参数根据被约束到离散图像网格位置上的特征来计算。如较早提及的那样，一些特征对没有被放置在精确对应的空间位置上。因此，所得到的集合可以包含逸出值和在反方向上使变换发生偏差的误差。在下述步骤中，Θ和C在子像素准确的迭代过程中被细化，以便实现更准确的对准。

利用J C和n[M来优化联合对应集合 $J=\{c_{i_1,j_1},c_{i_2,j_2},...,c_{i_n,j_n}\}$ 是值得期待的，该联合对应集合包含子像素准确对准的特征对并且在理想情况下不包含逸出值。优化后的联合对应集合的元素被用作ICP算法的输入，以便计算使全局图像相似性最大的变换：

$\hat{J}=\underset{J}{\arg \max }{L}_{\mathrm{global}}\left(J\right)=\underset{J}{\arg \max }\mathrm{ECC}(I_r,I_j^{T_i}).$

为了将特征对的数目保持在低值并保持配准效率，使用具有有限迭代步骤数的期望最大化(EM)类算法。在每次迭代中，从被逐步细化的联合对应集合J_k_J中计算出变换T_Jk。L_global被用作细化过程的收敛标准。一旦区域特征对已经以子像素精度被局部配准，这个特定对的下述配准并没有增强这个对应的质量，并且能够被忽略。因此，在每一个迭代步骤期间，通过只细化被迭代添加的特征对位置，计算时间能够被节省下来。

在图3中示出优化联合对应集合和配准变换的EM类算法的流程图。所述算法利用包含C中最前面两个对的联合对应集合： $J_0=\{c_{i_1,j_1},c_{i_2,j_2}\}$ 来初始化。对于这些初始对应，通常使用已排序集合的特征对中最好的两对，所述特征对在前一步中被获取。局部刚性配准被用来以子像素精度细化对应的显著区域特征。现在参考该图，在步骤31，子像素细化特征对应的当前集合被提供。在步骤32(估计步骤)，针对所有c_i，j∈C∧c_i，jJ^*，计算全局相似性度量L_global(J^*∪c_i，j)。在步骤33(最大化步骤)中，选择使L_global(J^*∪c_i，j)最大化的

然后，在步骤34，如果获得最大的L_global(J^*∪c_i，j)[L_global(J^*)，那么变换T_J*在步骤35被返回并且结束该算法。否则，在步骤36，利用局部刚性变换以子像素精度来对最大化特征对

进行配准，以细化特征中心。在步骤37，细化后的特征对

被增加到集合J^*： $J^{*}\&LeftArrow;J^{*}\&cup;c_{\hat{i},\hat{j}}^{*}$ 中，并且变换T_J*在步骤38中被重算。然后，步骤32-38被重复直到收敛。

根据本发明实施方式的算法已经对各种同一患者的三维医学图像进行了测试。已经对在不同时间获取的11个PET-CT体积对、3个不同治疗阶段的CT体积以及10个来自混合扫描仪的SPECT-CT体积对执行了测量。根据本发明实施方式的方法必须与不同模态、噪声、变化视野和一些PET-CT对中的图像强度伪像相竞争，其中一些切片具有在输入期间没有被校正的不同强度尺度。图4(a)-(c)是显示全部测量距离和标准偏差的表格，这些测量距离分别是针对以cm为单位给出的x、y和z方向上的PET-CT、CT-CT和SPECT-CT体积对。通过测量若干感兴趣点之间的距离由医学专家来评定PET-CT和CT-CT的配准质量，这些感兴趣点是：肺右尖和左尖，心尖，肝圆形边缘(liverr oundend)，左上和右上以及左下和右下肾缘。当10个SPECT-CT图像已经由现有技术的混合扫描仪获取时，医生严格随着x、y和z方向上从10到50mm的变化以及围绕每个轴从5到60度范围内变化的旋转手动撤销配准(deregister)SPECT图像。在配准之后，若干可识别的标志点已经在CT和SPECT图像上由医学专家选出。

已经在包含噪声或伪像的实时医学图像上进行了实验，该噪声或伪像是由于切片之间强度定标的变化。这些问题在配准之前没有被解决，以便用这样的数据来测试根据本发明实施方式的算法。图5(a)-(c)示出来自融合配准结果图像的三个切片，这些融合配准结果图像使用根据本发明的实施方式的算法获取到，并且分别来自PET-CT图像对、具有强度伪像的CT-CT图像对和SPECT-CT图像对。尽管在上述图像对中作为后者的CT图像使用有限视野来获取并且包含大量噪声，但是所提出的配准得到了可接受的准确度。剩余的误配准可以用本发明的非刚性变换的实施方式来解决。

医学专家使用用于可视化和测量的专用可视化软件来评定这些结果。针对该评估，医学专家可以在使用三维感兴趣区域的质心和标志点位置的直接标志点之间做出选择。通过将融合可视化和一些附加测量工具集成到再现软件中，这项任务得到了支持。在PET-CT情况下，z方向上的较高的标准偏差是显而易见的。出现这种情况的原因可能来源于采集模型之间的差别。CT图像显示出一次呼吸快照，而PET图像则在多个呼吸周期上被采集，并且或多或少地描述了平均呼吸运动。由于隔膜的这种运动，腹部区域中的一些器官被抬升和降低，这就导致了在数据采样中看到的较大的偏差。针对这个实验所使用的根据本发明实施方式的算法仅进行对刚性变换进行建模，而并没有对这样的局部变形进行建模。对于CT-CT数据，这种效应不再占优势，因为在两种采集中理想情况下患者吸气是相同的。SPECT-CT数据固有地匹配很好，并且SPECT上的用户定义的刚性变换并没有引入局部变形。因此，可以预期到对于这些情况会有良好的配准结果。

在所有结果中，因为医学专家必须通过点击各种切片视图中的位置来手动地指定位置，所以特定的测量误差被引入。但是，在已进行的这种类型的评估实验中，(观察者之间和同一观察者)在若干测量步骤中指定感兴趣点的距离的平均差不超过3mm。

应当理解，本发明能够用硬件、软件、固件、专用过程或者其组合的各种形式来实现。在一个实施方式中，本发明能够用软件实现为应用程序，该应用程序确实被包含在计算机可读的程序存储装置上。应用程序能够被上传到包括任何恰当结构的机器，并被该机器执行。

图7是用于实现根据本发明的实施方式的三维配准过程的示例性计算机***的框图。现在参考图7，用于实现本发明的计算机***71尤其可包括：中央处理单元(CPU)72，存储器73和输入/输出(I/O)接口74。计算机***71通常通过I/O接口74与显示器75和诸如鼠标和键盘的各种输入装置76连接。所述支持电路可包括诸如高速缓冲存储器、电源、时钟电路和通信总线的电路。存储器73可包括随机存取存储器(RAM)、只读存储器(ROM)、磁盘驱动器、磁带驱动器等或其组合。本发明能够被实施为例行程序77，该例行程序77被存储在存储器73中并通过CPU 72执行，以处理来自信号源78的信号。因而，该计算机***71是通用计算机***，该通用计算机***在执行本发明的例行程序77时成为专用计算机***。

所述计算机***71也包括操作***和微指令代码。在此所描述的各种过程和功能也可以是所述微指令代码的部分或通过操作***执行的应用程序的部分(或其组合)。另外，可以将各种其他***装置连接到该计算机平台上，诸如附加的数据存储装置和打印装置被连接到该计算机平台。

还应当理解，因为附图中所描述的一些组成***部件和方法步骤可以用软件实施，所以***部件(或过程步骤)之间的实际连接可以不同，这取决于本发明被编程的方式。给出在此所提供的本发明的教导，相关领域的技术人员能够设想到本发明的这些和类似的实施方式或配置。

虽然已参考优选实施方式对本发明进行了详细说明，但本领域技术人员应当理解，在不背离如所附权利要求中所阐述的本发明的精神和范围的情况下可对其进行各种修改和替换。

Claims (27)

1.一种对准一对图像的方法，该方法包括如下步骤：

提供具有第一图像和第二图像的一对图像，其中所述图像包括多个对应于三维空间中的像素域的强度；

在第一图像和第二图像中均标识显著特征区域，其中每个区域都与空间尺度相关；

通过每个区域的中心点来表示特征区域；

根据局部强度将一幅图像的特征点与另一幅图像的特征点进行配准；

通过类似性度量来对所述特征对进行排序；以及

通过将中心点细化到子像素精度来优化特征对的联合对应集合。

2.根据权利要求1所述的方法，还包括：在kD树中表示一幅图像的显著特征区域中心点；查询所述kD树的每个特征来找到一组最邻近的特征；以及从所述树中移除具有较低显著性值的那些最邻近特征和在所述每个特征的尺度内具有中心点的那些最邻近特征，其中获得显著特征区域在所述图像中的基本上均匀的分布。

3.根据权利要求1所述的方法，其中，所述空间尺度是包括所述特征区域的球形的半径。

4.根据权利要求2所述的方法，其中，所述kD树使用所述显著特征区域中心点的图像像素索引作为树叶，并且其中，从特征区域到最邻近特征区域的距离以图像索引为单位。

5.根据权利要求1所述的方法，其中，根据局部强度对特征点进行配准还包括：

使用所述第一图像和所述第二图像之间的迭代最近点变换来估计初始配准：

将所述第二图像的所有特征变换到所述第一图像的坐标空间中；

将所述所变换的特征存储在kD树中，并且根据预定的选择标准查询所述树的所述第一图像中的每个特征，以选择所述第二图像中的那些最邻近特征；和

在所述第一图像和所述第二图像中测试所述特征的所选择的特征对的平移不变性、旋转不变性和全局图像相似性度量，其中，通过所述所选择的特征对的全局图像相似性度量值来对所述所选择的特征对进行排序。

6.根据权利要求5所述的方法，其中，所述迭代最近点变换最小化了每组特征点之间的均方误差。

7.根据权利要求5所述的方法，其中，测试所述平移不变性包括估计 ${\widehat{\mathrm{\&Theta;}}}_{i,j}^T=p_i-p_j,$ 其中p_i和p_j是第i个第一图像特征和第j个第二图像特征在物理空间中的中心位置坐标。

8.根据权利要求5所述的方法，其中，测试所述旋转不变性包括估计

${\widehat{\mathrm{\&Theta;}}}_{i,j}^R=\underset{{\mathrm{\&Theta;}}^R}{\arg \max }\mathrm{ECC}(f_i,f_i^{T_{{\mathrm{\&Theta;}}^R}}),$

其中(f_i，f_j)分别表示所述第一图像和所述第二图像中的所述特征对，

$\mathrm{ECC}(f_i,f_j)=2-\frac{2H(f_i,f_j)}{H\left(f_i\right)+H\left(f_j\right)},$

其中H表示相对于在体素位置x和所述空间尺度s周围的球形邻近特征区域f_s内的图像强度值的熵，H被定义为

$H_D(s,x)=-{\&Integral;}_{R_s}p(i,s,x){\log }_{2}p(i,s,x)\mathrm{di},$

其中p(i，s，x)是被包含在f中的图像强度值i的概率密度函数，其中H(f_i，f_j)是联合微分熵，H(f_i，f_j)被定义为：

$H(f_i,f_j)=-{\&Integral;}_{f_i,f_j}p(f_i,f_j){\log }_{2}p(f_i,f_j)\mathrm{dIdJ},$

其中p(f_i，f_j)是特征区域f_i和f_j中的图像强度的联合概率密度，并且I和J分别在所述第一和第二图像中的可能的强度值的集合中取值。

9.根据权利要求5所述的方法，其中，测试所述全局图像相似性度量包括估计

$L_{\mathrm{global}}\left(c_{i,j}\right)=\mathrm{ECC}(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}),$

其中I_r表示所述第一图像，

表示所述第二图像到所述第一图像的坐标空间上的所述变换，

$\mathrm{ECC}(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})=2-\frac{2H(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})}{H\left(I_r\right)+H\left(H_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}\right)},$

其中H表示相对于在体素位置x和所述空间尺度s周围的所述图像中的一幅图像内的图像强度值的熵，H被定义为：

H(s，x)＝-∫_Ip(i，s，x)log₂p(i，s，x)di，

其中p(i，s，x)是被包含在I中的图像强度值i的概率密度函数，其中是联合微分熵，

被定义为：

$H(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})=-{\&Integral;}_{I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}}p(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}){\log }_{2}p(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})\mathrm{dIdJ},$

其中

是图像I_r和

中的图像强度的联合概率密度，并且I和J分别在所述第一和第二图像中的可能的强度值的集合中取值，并且其中L_global在所述第一和第二图像的整个重叠域上被评估。

10.根据权利要求1所述的方法，其中，优化特征对的联合对应集合还包括：

利用根据所述相似性度量为最相似的特征对来对所述联合对应集合进行初始化；

针对所述联合对应集合与没有已被包括在所述联合对应集合中的每个特征对的并集估计所述相似性度量；

选择最大化所述并集的相似性度量的特征对，其中，如果所述最大化特征对与所述联合对应集合的并集的相似性度量大于该联合对应集合的相似性度量，则所述最大化特征对利用局部刚性变换以子像素精度进行配准并且被增加到所述联合对应集合中。

11.根据权利要求10所述的方法，其中，通过使用迭代最近点过程计算特征对之间的配准变换来使所述相似性度量最大化。

12.根据权利要求10所述的方法，其中，如果所述最大化特征对与所述联合对应集合的并集的相似性度量小于或等于该联合对应集合的相似性度量，则提供配准变换，所述配准变换由最大化所述相似性度量的特征对之间的配准变换计算得到。

13.一种对准一对图像的方法，该方法包括如下步骤：

提供具有第一图像和第二图像的一对图像，其中，所述图像包括多个对应于三维空间中的像素域的强度；

在第一图像和第二图像中均标识显著特征区域，其中每个区域都与空间尺度相关；

使用所述第一图像和所述第二图像之间的迭代最近点变换来估计初始配准；

将所述第二图像的所有特征变换到所述第一图像的坐标空间；

将所述所变换的特征存储在kD树中，并且根据预定的选择标准来查询所述树的所述第一图像中的每个特征，以选择所述第二图像中的那些最邻近特征；

在所述第一图像和所述第二图像中测试所述特征的所选择的特征对的平移不变性、旋转不变性和全局图像相似性度量；以及

通过所述所选择的特征对的全局图像相似性度量值来对所述所选择的特征对进行排序。

14.根据权利要求13所述的方法，还包括：通过每个区域的中心点表示特征区域；将一幅图像的所述特征区域中心点存储在kD树中；查询所述kD树的每个特征来找到一组最邻近特征；以及从所述树中移除具有较低显著性值的那些最邻近特征和在所述每个特征的尺度内具有中心点的那些最邻近特征，其中，获得显著特征区域在所述图像中的基本上均匀的分布。

15.根据权利要求13所述的方法，还包括：

利用根据所述相似性度量为最相似的特征对来对联合对应集合进行初始化；

针对所述联合对应集合与没有已被包括在所述联合对应集合中的每个特征对的并集估计所述相似性度量；

选择最大化所述并集的相似性度量的特征对，其中，如果所述最大化特征对与所述联合对应集合的并集的相似性度量大于该联合对应集合的相似性度量，则所述最大化特征对利用局部刚性变换以子像素精度进行配准并且被增加到所述联合对应集合中；

其中，所述全局图像相似性度量被定义为

$L_{\mathrm{global}}\left(c_{i,j}\right)=\mathrm{ECC}(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}),$

其中I_r表示所述第一图像，

表示所述第二图像到所述第一图像的坐标空间上的所述变换，

$\mathrm{ECC}(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})=2-\frac{2H(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})}{H\left(I_r\right)+H\left(H_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}\right)},$

其中H表示相对于在体素位置x和所述空间尺度s周围的所述图像中的一幅图像内的图像强度值的熵，H被定义为：

H(s，x)＝-∫_Ip(i，s，x)log₂p(i，s，x)di，

其中p(i，s，x)是被包含在I中的图像强度值i的概率密度函数，其中是联合微分熵，

被定义为：

$H(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})=-{\&Integral;}_{I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}}p(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}){\log }_{2}p(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})\mathrm{dIdJ},$

其中

是图像I_r和

中的图像强度的联合概率密度，并且I和J分别在所述第一和第二图像中的可能的强度值的集合中取值，并且其中L_global在所述第一和第二图像的整个重叠域上被评估。

16.一种计算机可读程序存储装置，其确实地包含了可由计算机执行的程序指令，以执行用于对准一对图像的方法步骤，所述方法包括如下步骤：

提供具有第一图像和第二图像的一对图像，其中所述图像包括多个对应于三维空间中的像素域的强度；

在第一图像和第二图像中均标识显著特征区域，其中每个区域都与空间尺度相关；

通过每个区域的中心点来代表特征区域；

根据局部强度将一幅图像的特征点与另一幅图像的特征点进行配准；

通过类似性度量来对所述特征对进行排序；以及

通过将中心点细化到子像素精度来优化特征对的联合对应集合。

17.根据权利要求16所述的计算机可读程序存储装置，所述方法还包括：在kD树中表示一幅图像的显著特征区域中心点；查询所述kD树的每个特征来找到一组最邻近特征；以及从所述树中移除具有较低显著性值的那些最邻近特征和在所述每个特征的尺度内具有中心点的那些最邻近特征，其中，获得显著特征区域在所述图像中的基本上均匀的分布。

18.根据权利要求16所述的计算机可读程序存储装置，其中，所述空间尺度是包括所述特征区域的球形的半径。

19.根据权利要求17所述的计算机可读程序存储装置，其中，所述kD树使用所述显著特征区域中心点的图像像素索引作为树叶，并且其中，从特征区域到最邻近特征区域的距离以图像索引为单位。

20.根据权利要求16所述的计算机可读程序存储装置，其中，根据局部强度对特征点进行配准还包括：

使用所述第一图像和所述第二图像之间的迭代最近点变换来估计初始配准：

将所述第二图像的所有特征变换到所述第一图像的坐标空间中；

将所述所变换的特征存储在kD树中，并且根据预定的选择标准查询所述树的所述第一图像中的每个特征，以选择所述第二图像中的那些最邻近特征；和

在所述第一图像和所述第二图像中测试所述特征的所选择的特征对的平移不变性、旋转不变性和全局图像相似性度量，其中，通过所述所选择的特征对的全局图像相似性度量值来对所述所选择的特征对进行排序。

21.根据权利要求20所述的计算机可读程序存储装置，其中，所述迭代最近点变换最小化每组特征点之间的均方误差。

22.根据权利要求20所述的计算机可读程序存储装置，其中，测试所述平移不变性包括估计 ${\widehat{\mathrm{\&Theta;}}}_{i,j}^T=p_i-p_j,$ 其中p_i和p_j是第i个第一图像特征和第j个第二图像特征在物理空间中的中心位置坐标。

23.根据权利要求20所述的计算机可读程序存储装置，其中，测试所述旋转不变性包括估计

${\widehat{\mathrm{\&Theta;}}}_{i,j}^R=\underset{{\mathrm{\&Theta;}}^R}{\arg \max }\mathrm{ECC}(f_i,f_i^{T_{{\mathrm{\&Theta;}}^R}}),$

其中(f_i，f_j)分别表示所述第一图像和所述第二图像中的所述特征对，

$\mathrm{ECC}(f_i,f_j)=2-\frac{2H(f_i,f_j)}{H\left(f_i\right)+H\left(f_j\right)},$

其中H表示相对于在体素位置x和所述空间尺度s周围的球形邻近特征区域f_s内的图像强度值的熵，H被定义为

$H_D(s,x)=-{\&Integral;}_{R_s}p(i,s,x){\log }_{2}p(i,s,x)\mathrm{di},$

其中p(i，s，x)是被包含在f中的图像强度值i的概率密度函数，其中H(f_i，f_j)是联合微分熵，H(f_i，f_j)被定义为：

$H(f_i,f_j)=-{\&Integral;}_{f_i,f_j}p(f_i,f_j){\log }_{2}p(f_i,f_j)\mathrm{dIdJ},$

其中p(f_i，f_j)是特征区域f_i和f_j中的图像强度的联合概率密度，并且I和J分别在所述第一和第二图像中的可能的强度值的集合中取值。

24.根据权利要求20所述的计算机可读程序存储装置，其中，测试所述全局图像相似性度量包括估计

$L_{\mathrm{global}}\left(c_{i,j}\right)=\mathrm{ECC}(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}),$

其中I_r表示所述第一图像， 表示所述第二图像到所述第一图像的坐标空间上的所述变换，

$\mathrm{ECC}(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})=2-\frac{2H(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})}{H\left(I_r\right)+H\left(H_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}\right)},$

其中H表示相对于在体素位置x和所述空间尺度s周围的所述图像中的一幅图像内的图像强度值的熵，其被定义为：

H(s，x)＝-∫_Ip(i，s，x)log₂p(i，s，x)di，

其中p(i，s，x)是被包含在I中的图像强度值i的概率密度函数，其中是联合微分熵， 被定义为：

$H(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})=-{\&Integral;}_{I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}}p(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}}){\log }_{2}p(I_r,I_t^{T_{{\widehat{\mathrm{\&Theta;}}}_{i,j}}})\mathrm{dIdJ},$

其中 是图像I_r和 中的图像强度的联合概率密度，并且I和J分别在所述第一和第二图像中的可能的强度值的集合中取值，并且其中L_global在所述第一和第二图像的整个重叠域上被评估。

25.根据权利要求16所述的计算机可读程序存储装置，其中，优化特征对的联合对应集合还包括：

利用根据所述相似性度量为最相似的特征对来对所述联合对应集合进行初始化；

针对所述联合对应集合与没有已被包括在所述联合对应集合中的每个特征对的并集估计所述相似性度量；

选择最大化所述并集的相似性度量的特征对，其中，如果所述最大化特征对与所述联合对应集合的并集的相似性度量大于该联合对应集合的相似性度量，则所述最大化特征对利用局部刚性变换以子像素精度进行配准并且被增加到所述联合对应集合中。

26.根据权利要求25所述的计算机可读程序存储装置，其中，通过使用迭代最近点过程计算特征对之间的配准变换来使所述相似性度量最大化。

27.根据权利要求25所述的计算机可读程序存储装置，其中，如果所述最大化特征对与所述联合对应集合的并集的相似性度量小于或等于该联合对应集合的相似性度量，则提供配准变换，所述配准变换由最大化所述相似性度量的特征对之间的配准变换计算得到。

Images