CN106127726B - 一种非特征提取和非参数的3d图像配准方法 - Google Patents

Abstract

本发明涉及图像融合和图像序列监测等技术领域，具体涉及了一种基于非特征提取和非参数3D图像配准方法。该方法用边缘检测器来检测要观察的参考图像Z_R的边缘体素，对给定参考图像R的体素(x,y,z)进行判断，对R局部无退化体进行几何变换，用R的局部退化点定义T(x,y,z)，之后输出几何变换参数。该方法由于局部平滑的灵活性，该方法相关的几何变换不需要任何参数形式，并且在实践中实用有效。

Description

一种非特征提取和非参数的3D图像配准方法

技术领域

本发明涉及图像融合和图像序列监测等技术领域，具体涉及了一种基于非特征提取和非参数3D图像配准方法。

背景技术

在许多图像应用中，为了完成图片对象的融合和差异性的检测，我们都需要比对同一个对象的多张图片。比如，为了密切监测脑肿瘤患者的肿瘤生长情况，我们需要对该患者在不同的时刻拍的脑部的核磁共振图像进行比对。因为采集图像时，成像设备和图像对象的相对位置在不同的时间几乎不会相同。为了能让比对有效，我们需要快速地完成相关图像的几何匹配。图像配准(IR)就是为了达到这个目的,它是包括医疗成像、遥感技术、指纹和脸部识别、图像压缩、视频增强等许多图像应用在内的基础任务。

文献中存在的大部分的图像配准(IR)方法都是对2D图像的分析，他们大致可以分为基于特征和基于灰度方法的两类。要使用基于特征的方法完成两张图像的匹配，我们首先要考虑选取两张图像的两个特征集，然后去寻找能够很好的匹配两个特征集的几何变换。特征提取通常是一个耗时多随机性强且具有挑战性的任务。

因为图像获取技术的快速发展，3D图像在某些应用中越来越受欢迎，比如医疗成像。因此，随着3D图像的广泛引用，3D图像配准成为一项很重要的研究。3D图像的结构要比2D图像的复杂得多，例如，3D图像的边缘位置的表面结构比2D图像的边缘曲线复杂。到目前为止，在文献中也没有太多的3D图像配准的方法。基于迭代最近点(ICP)方法是种常用的3D配准方法，这个基于特征的IR方法是建立在有参数形式的几何变换的假设上。要用这个方法，首先要从两个相关图像提取两个特征集，然后ICP方法会在估计几何变换的参数和两张图片间最佳点对应交替进行。显然，这个方法有这几个缺点：(1)提取相关的特征集很不方便；(2)利用参数化几何变换的限制；(3)严重依赖方法的初始化质量；(4)容易获得局部极小，而不是全局极小。

发明内容

本发明提出了一个新的3D IR方法，这个方法用在一个图像到另一个图像的几何变换是非刚性体的情况中。刚体变换意味着任何两个体素在图像之间的欧氏距离在几何变换后不会改变。而非刚体变换更灵活，它可以应对关节物体或软物体随着时间的推移改变形状的情况，典型的非刚体IR方法常应用在软组织(例如，人体器官)等很常见的生物医学研究中。在提出的新方法中，我们采用非参数局部平滑方法，通过这样的方法，在一个给定的体素的估计仅取决于3D图像中在它周围的图像灰度。这个局部的估计性质使得它可以在几何变换中不强加限制性的假设。

3D图像配准问题描述如下：设R和M为两张要配准的3D图像，R通常被称为参考图像，M为移动图像。他们真实的灰度函数可表示为R(x,y,z)和M(x,y,z)，假设他们有如下的关系：

M(T₁(x,y,z),T₂(x,y,z),T₃(x,y,z))＝R(x,y,z),(x,y,z)∈Ω (1)

Ω是图像R的设计空间，T(x,y,z)＝(T₁(x,y,z),T₂(x,y,z),T₃(x,y,z))是未知的几何变换估计，设计的方法通过检测到的两张图像的图像灰度的统计模型估计T(x,y,z)。

Z_M(x_i,y_j,z_k)＝M(x_i,y_j,z_k)+εM(x_i,y_j,z_k)

Z_R(x_i,y_j,z_k)＝R(x_i,y_j,z_k)+εR(x_i,y_j,z_k)

i,j,k＝1,2,....,n (2)

{(x_i,y_j,z_k)}是体素的位置，εM(x_i,y_j,z_k)和εR(x_i,y_j,z_k)分别是两幅图像中均值为0、方差未知、独立同分布的随机误差。非参数的IR方法通过检测的图像灰度去估计T(x,y,z)，并不需要对T(x,y,z)加上参数形式。现给定一个三维体素(x,y,z)∈Ω，那么我们可以写出：

其中，b(x,y,z)＝T₁(x,y,z)-x,c(x,y,z)＝T₂(x,y,z)-y，d(x,y,z)＝T₃(x,y,z)-z。因此，对T(x,y,z)就等于对(b(x,y,z),c(x,y,z),d(x,y,z))的估计。估计值就可以写为：

现给定三维体素(x,y,z)∈Ω，如果几何变换的幅度T(x,y,z)很小，并且M具有一级偏导数在(x,y,z)，由泰勒的展开式，有：

M(T₁(x,y,z),T₂(x,y,z),T₃(x,y,z))＝M(x,y,z)+M'_x(x,y,z)b(x,y,z)

+M'_y(x,y,z)c(x,y,z)+M'_z(x,y,z)d(x,y,z)+ο(||T(x,y,z)-(x,y,z)||)

(4)

||·||是欧几里得范数。在这种情况下，通过公式(1)和公式(4)，R(x,y,z)＝M(T₁(x,y,z),T₂(x,y,z),T₃(x,y,z))可以近似为M(x,y,z)+M'_x(x,y,z)b(x,y,z)+M'_y(x,y,z)c(x,y,z)+M'_z(x,y,z)d(x,y,z)。因此，(b(x,y,z),c(x,y,z),d(x,y,z))可以当成是逼近误差。

R(x,y,z)-[M(x,y,z)+M'_x(x,y,z)b(x,y,z)+M'_y(x,y,z)c(x,y,z)+M'_z(x,y,z)d(x,y,z)]要尽可能的小。实际中，R(x,y,z)，M(x,y,z)，M'_x(x,y,z)，M'_y(x,y,z)，M'_z(x,y,z)都是难以观察的。所观察到的图像强度{Z_M(x_i,y_j,z_k)}和{Z_R(x_i,y_j,z_k)}定义在含有随机噪声的公式(2)。为了能在平滑噪声的同时也能估计出(b(x,y,z),c(x,y,z),d(x,y,z))，我们在下面采用局部线性核(LLK)估计的思想统计非参数回归。首先，M'_x(x,y,z)，M'_y(x,y,z)，M'_z(x,y,z)的值可以分别通过他们传统的LLK估计值 估计出来。然后，在一个半径为h_n的球形邻域(x,y,z)可表示为Ο(x,y,z；h_n)，b(x,y,z)，c(x,y,z)和d(x,y,z)可以通过以下的最小化问题的解决方案来估计：

是一个由单位圆支撑的三元密度核函数。求解b(x,y,z),c(x,y,z)和d(x,y,z)最小估计值要能满足逼近误差平方和最小和内核函数K确定的权重。在统计文献中，K经常被选为Epanechnikov核函数C_K(1-x²)(1-y²)(1-z²)I(x²+y²+z²≤1)，C_K是归一化常数，I(·)是指示函数。公式(5)中，如果一个三维体素(x_i,y_j,z_k)∈Ο(x,y,z；h_n)和给定的三维体素(x,y,z)距离很远，那么(x_i,y_j,z_k)对应的逼近误差会有一个很小的权重。

这也就不难去验证公式(6)给出了问题(5)给出的解决方法。公式6中，当s₁,s₂＝x,y,z

M'_x(x,y,z)，M'_y(x,y,z)，M'_z(x,y,z)传统的LLK估计值为：

当s＝x,y,z

从上面的描述中，我们知道公式(6)是在如下理想条件中获得：(1)||T(x,y,z)-(x,y,z)||足够小，能够满足公式(4)中M(T₁(x,y,z),T₂(x,y,z),T₃(x,y,z))良好的一阶近似。(2)M在(x,y,z)中有一阶偏导数。(3)公式(6)中等式右边的分母部分不为零。上面的条件(1)和条件(2)表明，公式(3)和公式(6)定义的估计值在变换比较大或者M不平滑(例如，M的边缘位置)的情况下，不能很好地估计T(x,y,z)。条件(3)则表明，在如下的等式中，估计值不能被很好地定义。

从数学上来讲，它能够证明：(1)当M具有连续的一阶衍生性时，随着n值的增大，会分别逼近M'_x(x,y,z)，M'_y(x,y,z)，M'_z(x,y,z)。(2)如果在公式(7)中的 分别被M'_x(x,y,z)，M'_y(x,y,z)，M'_z(x,y,z)替换，那么M在邻域Ο(x,y,z；h_n)中当且仅当有如下的一个单变量连续可微函数和一个常量ρ才能满足公式(7)的等式。

对所有的(x',y',z')∈Ο(x,y,z；h_n)都成立。

如果M在Ο(x,y,z；h_n)中满足公式(8)，它的灰度等级和Ο(x,y,z；h_n)中的线段ρx'+y'＝ρ₀是相同的，所有能够满足这条线段的常量ρ₀都包含在Ο(x,y,z；h_n)中。在这种情况下，三元函数M就能在Ο(x,y,z；h_n)简化为局部性的，并且它也不可能唯一的确定T(x,y,z)，因为当(x',y',z')∈Ο(x,y,z；h_n)时，任何沿线方向的微小移动都能改变M(x',y',z')的值。在本文中，如果一个三维体素(x,y,z)∈Ω，M在这个三维体素(x,y,z)有偏导数，并且存在能够满足公式(8)的邻域Ο(x,y,z；h_n)，那么这样一个三维体素被称为局部退化的三维体素。同样的，我们可以给参考图像R定义局部退化体素和局部不退化体素。

只有在局部无退化体素中或者当参考图像R的图像灰度函数不平滑时，几何变换T(x,y,z)才能被正确的估计。基于这样的原因，我们提出的方法包括五个主要的步骤，具体介绍如下：

步骤1：用边缘检测器来检测要观察的参考图像Z_R的边缘体素；

步骤2：对于已经给定的真实的参考图像R的体素(x,y,z)，以该体素为中心，半径为r₁的圆形邻域表示为Ο(x,y,z；r₁)，如果在Ο(x,y,z；r₁)中检测到的边缘体素数量小于[nr₁]，其中，[]为取整运算符号，那么(x,y,z)就被认为是参考图像R的连续体素，在这种情况下，如果公式(6)等式右边的分母(在M被R替换之后)大于或者等于预先指定的阈值μ_n，那么(x,y,z)就被认为是R的局部无退化体素；否则，它就是局部退化体素；

步骤3：设D为R局部无退化体素或在半径为r₁邻域内至少含有[nr₁]边缘体素的体素的集合，对于所有的(x,y,z)∈D，几何变换的幅度T(x,y,z)都将通过下面的步骤(3)和步骤(4)的过程中方法进行计算，对于移动图像上的任意体素(x',y',z')∈O(x,y,z；r₁)，只考虑它的圆形邻域Ο(x',y',z'；r₂)，这个邻域的半径r₂可以和r₁相同，也可以不同，计算平均平方差MSD有：

是Ο(x',y',z'；r₂)的体素数量，被定义为下面最小化的表示：

步骤4：如果(x,y,z)是R的局部退化点，那么T(x,y,z)的估计值被定义如下：首先，在D中找到最接近(x,y,z)的体素，表示为(x⁽¹⁾,y⁽¹⁾,z⁽¹⁾)；然后，令

如果则否则定义，

步骤5：输出几何变换参数。

本发明的有益效果是：该方法对初始化质量的依赖性降低，该方法相关的几何变换不需要任何参数形式；实践中更有效。

附图说明

图1：一种非特征提取和非参数信息的图像配准方法流程图。

具体实施方式

下面结合附图和实例对发明作进一步说明：

我们提出的图像配准程序包括五个主要的步骤，具体介绍如下：

步骤1：用边缘检测器来检测要观察的参考图像Z_R的边缘体素；

步骤2：对于已经给定的真实的参考图像R的体素(x,y,z)，以该体素为中心，半径为r₁的圆形邻域表示为Ο(x,y,z；r₁)，如果在Ο(x,y,z；r₁)中检测到的边缘体素数量小于[nr₁]，其中，[]为取整运算符号，那么(x,y,z)就被认为是参考图像R的连续体素，在这种情况下，如果公式(6)等式右边的分母(在M被R替换之后)大于或者等于预先指定的阈值μ_n，那么(x,y,z)就被认为是R的局部无退化体素；否则，它就是局部退化体素；

步骤3：设D为R局部无退化体素或在半径为r₁邻域内至少含有[nr₁]边缘体素的体素的集合，对于所有的(x,y,z)∈D，几何变换的幅度T(x,y,z)都将通过下面的步骤(3)和步骤(4)的过程中方法进行计算，对于移动图像上的任意体素x',y',z'∈O(x,y,z；r₁)，只考虑它的圆形邻域Ο(x',y',z'；r₂)，这个邻域的半径r₂可以和r₁相同，也可以不同，计算平均平方差MSD有：

是Ο(x',y',z'；r₂)的体素数量，被定义为下面最小化的表示：

步骤4：如果(x,y,z)是R的局部退化点，那么T(x,y,z)的估计值被定义如下：首先，在D中找到最接近(x,y,z)的体素，表示为(x⁽¹⁾,y⁽¹⁾,z⁽¹⁾)；然后，令

如果则否则定义，

步骤5：输出几何变换参数。

Claims (2)

1.一种基于非特征提取和非参数3D图像配准方法，其特征在于包括以下步骤：

(1)用边缘检测器来检测要观察的参考图像Z_R的边缘体素；

(2)对于已经给定的真实的参考图像R的体素(x,y,z)，以该体素为中心，半径为r₁的圆形邻域表示为Ο(x,y,z；r₁)，如果在Ο(x,y,z；r₁)中检测到的边缘体素数量小于[nr₁]，其中，[]为取整运算符号，那么(x,y,z)就被认为是参考图像R的连续体素，在这种情况下，继续判断体素是局部无退化体素还是局部退化体素；

(3)设D为R局部无退化体素或在半径为r₁邻域内至少含有[nr₁]边缘体素的体素的集合，对于所有的(x,y,z)∈D，几何变换的幅度T(x,y,z)都将通过下面的步骤(3)和步骤(4)的过程中方法进行计算，对于移动图像上的任意体素(x',y',z')∈O(x,y,z；r₁)，只考虑它的圆形邻域Ο(x',y',z'；r₂)，这个邻域的半径r₂可以和r₁相同，也可以不同，计算平均平方差MSD有：

是Ο(x',y',z'；r₂)的体素数量，是让下式最小的几何变换：

(4)如果(x,y,z)是R的局部退化点，那么T(x,y,z)的估计值被定义如下：首先，在D中找到最接近(x,y,z)的体素，表示为(x⁽¹⁾,y⁽¹⁾,z⁽¹⁾)；然后，令

如果则否则定义，

(5)输出几何变换参数。

2.根据权利要求1所述的一种基于非特征提取和非参数3D图像配准方法，其特征在于：在权利要求1的步骤(3)中，MSD被用作匹配准则，或者使用包括交叉准则和图像差熵在内的准则作为匹配准则。