CN107123133A

CN107123133A - 非特征的3d图像快速刚性配准方法

Info

Publication number: CN107123133A
Application number: CN201710300743.4A
Authority: CN
Inventors: 宋华军; 任鹏; 吴玉兴; 于玮; 钟鑫
Original assignee: China University of Petroleum East China
Current assignee: China University of Petroleum East China
Priority date: 2017-05-02
Filing date: 2017-05-02
Publication date: 2017-09-01

Abstract

本发明提供了一种非特征的3D图像快速刚性配准方法。该非特征的3D图像快速刚性配准方法，采用泰勒展开式和最小二乘法，通过迭代计算获得收敛时的几何变换参数，实现了3D图像快速刚性配准，能够获得较高的精度和速度，可以在3D图像采集数据大幅减少的情况下进行有效计算。

Description

非特征的3D图像快速刚性配准方法

技术领域

本发明涉及图像融合和图像序列监测技术领域，具体涉及一种非特征的3D图像快速刚性配准方法。

背景技术

3D图像在应用中越来越受欢迎。例如在MRI(磁共振成像)中，人们习惯于从3D物体中获取一组3D图像来研究3D物体(如：患者的头部)的生物机制。为了合理地去比较同一个对象在不同时刻两个3D图像的差异，需要事先对其进行配准(IR)。图像配准技术广泛的应用在医疗成像，目标跟踪，遥感技术，指纹和脸部识别，图像压缩和视频增强等场合。

文献中存在的大部分的2D图像配准方法大致可以分为基于特征和基于灰度方法的两类。基于特征的方法首先要考虑选取两张图像的两个特征集，然后去寻找能够很好的匹配两个特征集的几何变换F。常用的特征包括角点、LBP、SURF、质心或模板，还有图像灰度函数得到的退化像素。由于特征提取通常是一个耗时多、随意性强且具有挑战性的任务，最近IR研究更多关注直接对两幅图像灰度值的变换F，包括基于参数的转换，和更灵活的非参数的转换。

一些2D图像配准方法被引入到3D Slicer软件包(http://slicer.org/)中从而推广到3D情况下。最近比较流行的ICP算法是基于点云数据集，点云数据是3D扫描设备产生的，主要用来代表物体的外表面形状，而现实中很多的3D图像是存放整个图像的各个位置的灰度信息(比如MRI图像)，基于点云的配准算法不适用这些情况。

综上可知，现有的一些配准方法在用于3D图像配准时，往往存在配准结果不准确的问题。

发明内容

在下文中给出了关于本发明的简要概述，以便提供关于本发明的某些方面的基本理解。应当理解，这个概述并不是关于本发明的穷举性概述。它并不是意图确定本发明的关键或重要部分，也不是意图限定本发明的范围。其目的仅仅是以简化的形式给出某些概念，以此作为稍后论述的更详细描述的前序。

鉴于此，本发明提出一种非特征的3D图像快速刚性配准方法，以至少解决现有配准方法在用于3D图像配准时存在配准结果不准确的问题。

根据本发明的一个方面，提供了一种非特征的3D图像快速刚性配准方法，所述非特征的3D图像快速刚性配准方法包括：步骤一、获得实际参考图像Z_R(x_i,y_j,z_k)和实际移动图像Z_M(x_i,y_j,z_k)；执行步骤二；步骤二、估计下式中的和

其中，和分别表示移动图像f_M(x,y,z)在(x,y,z)上的一阶偏导数f′_Mx(x,y,z)、f′_My(x,y,z)和f′_Mz(x,y,z)的估计函数；(x_i,y_j,z_k)∈Ω表示参考图像与移动图像对应位置的像素点，i,j,k指的是各像素点的位置下标，Z_M(x_i,y_j,z_k)表示的是实际移动图像，K_h是以单位圆支撑的核函数，K_h(x,y,z)＝(1-x²)(1-y²)(1-z²)，h＞0是一个带宽参数；执行步骤三；步骤三、由Θ＝(X^TX)^-1X^TY获取参数Θ的当前估计值其中，Θ＝(α,β,γ,Δx,Δy,Δz)^T，其中，X是n³×6设计矩阵，它的第[(i-1)n²+(j-1)n+k]行元素是由确定的；Y是n³维向量，它的第[(i-1)n²+(j-1)n+k]元素由Z_R(x_i,y_j,z_k)-Z_M(x_i,y_j,z_k)确定；根据参数Θ的当前估计值确定几何变换函数F的当前估计值其中，执行步骤四；步骤四、判定当前是否满足以及若是，执行步骤六；否则，执行步骤五；其中，ε₁和ε₂为0至1之间的实数；步骤五、根据几何变换函数F的当前估计值对当前Z_M(x_i,y_j,z_k)进行变换，将Z_M(x_i,y_j,z_k)变换后的值记为并用更新Z_M(x_i,y_j,z_k)后返回执行步骤二；步骤六、将参数Θ的当前估计值确定为计算结果，以完成图像配准。

进一步地，ε₁＝ε₂＝0.1。

进一步地，在定义几何变换函数F(x,y,z)时，(x,y,z)是f_M(x,y,z)的像素点；||F(x,y,z)-(x,y,z)||小于预设实数，以使f_M(F₁(x,y,z),F₂(x,y,z),F₃(x,y,z))的泰勒展开式有效。

现有技术中，大都基于特征点和针对特定的约束条件，带来了特征选择耗时多，随机性强，而且约束条件使用不灵活等问题。相比之下，本发明的非特征的3D图像快速刚性配准方法针对这些问题，提出直接使用图像灰度值的无特征3D刚性配准方法，该方法使用泰勒展开式和最小二乘法直接计算待配准图像的变换参数，并且使用较少的数据点完成快速的配准。本发明的算法能够获得较高的精度，并且在数据减少时仍可以有效计算；而且，该方法在需要使用数据压缩的应用里是非常实用的，提出的方法适合于应用在当数据量较大或者需要快速进行图像配准场合。

通过以下结合附图对本发明的最佳实施例的详细说明，本发明的这些以及其他优点将更加明显。

附图说明

本发明可以通过参考下文中结合附图所给出的描述而得到更好的理解，其中在所有附图中使用了相同或相似的附图标记来表示相同或者相似的部件。所述附图连同下面的详细说明一起包含在本说明书中并且形成本说明书的一部分，而且用来进一步举例说明本发明的优选实施例和解释本发明的原理和优点。在附图中：

图1是示意性地示出本发明的非特征的3D图像快速刚性配准方法的一个示例的结构图。

本领域技术人员应当理解，附图中的元件仅仅是为了简单和清楚起见而示出的，而且不一定是按比例绘制的。例如，附图中某些元件的尺寸可能相对于其他元件放大了，以便有助于提高对本发明实施例的理解。

具体实施方式

在下文中将结合附图对本发明的示范性实施例进行描述。为了清楚和简明起见，在说明书中并未描述实际实施方式的所有特征。然而，应该了解，在开发任何这种实际实施例的过程中必须做出很多特定于实施方式的决定，以便实现开发人员的具体目标，例如，符合与系统及业务相关的那些限制条件，并且这些限制条件可能会随着实施方式的不同而有所改变。此外，还应该了解，虽然开发工作有可能是非常复杂和费时的，但对得益于本公开内容的本领域技术人员来说，这种开发工作仅仅是例行的任务。

在此，还需要说明的一点是，为了避免因不必要的细节而模糊了本发明，在附图中仅仅示出了与根据本发明的方案密切相关的装置结构和/或处理步骤，而省略了与本发明关系不大的其他细节。

本发明的实施例提供了一种非特征的3D图像快速刚性配准方法，所述非特征的3D图像快速刚性配准方法包括：步骤一、获得实际参考图像Z_R(x_i,y_j,z_k)和实际移动图像Z_M(x_i,y_j,z_k)；执行步骤二；步骤二、估计下式中的和

图1给出了本发明的非特征的3D图像快速刚性配准方法的一个示例性处理。

如图1所示，该方法开始后，首先执行步骤一。

在步骤一中，获得实际参考图像ZR(xi,y_j,zk)和实际移动图像ZM(xi,y_j,zk)，然后，执行步骤二。

在步骤二中，估计公式一一至公式一三中的和然后，执行步骤三。

公式一一：

公式一二：

公式一三：

其中，和分别表示移动图像f_M(x,y,z)在(x,y,z)上的一阶偏导数f′_Mx(x,y,z)、f′_My(x,y,z)和f′_Mz(x,y,z)的估计函数；(x_i,y_j,z_k)∈Ω表示参考图像与移动图像对应位置的像素点，i,j,k指的是各像素点的位置下标，Z_M(x_i,y_j,z_k)表示的是实际移动图像，K_h是以单位圆支撑的核函数，K_h(x,y,z)＝(1-x²)(1-y²)(1-z²)，h＞0是一个带宽参数。然后，执行步骤三。

在步骤三中，由公式二获取参数Θ的当前估计值

公式二：Θ＝(X^TX)^-1X^TY。

其中，Θ＝(α，β，γ，Δx，Δy，Δz)^T，Θ是几何变换函数F的参数(即几何变换参数)。

其中，在参数Θ中，α、β和γ分别表示沿x轴、y轴和z轴的旋转角度，而Δx、Δy和Δz分别表示沿着x轴、y轴和z轴这三个轴的平移量。

其中，X是n³×6设计矩阵，它的第[(i-1)n²+(j-1)n+k]行元素是由公式三确定的。

公式三：

此外，Y是n³维向量，它的第[(i-1)n²+(j-1)n+k]元素由Z_R(x_i,y_j,z_k)-Z_M(x_i,y_j,z_k)确定

这样，可以根据参数Θ的当前估计值确定几何变换函数F的当前估计值几何变换函数F可由公式四一至公式四四确定，相应地，几何变换函数F的估计值可由公式五一至公式五四确定。然后，执行步骤四。

公式四一：F(x,y,z)＝(F₁(x,y,z),F₂(x,y,z),F₃(x,y,z))

公式四二：F₁(x,y,z)＝x(cosαcosβ)+y(sinαcosβ)-z(sinβ)+Δx

公式四三：

F₂(x,y,z)＝x(cosαsinβsinγ-sinαcosγ)+y(sinαsinβsinγ+cosαcosγ)+z(cosβsinγ)+Δy

公式四四：

F₃(x,y,z)＝x(cosαsinβcosγ+sinαsinγ)+y(sinαsinβcosγ-cosαsinγ)+z(cosβcosγ)+Δz

公式五一：

公式五二：

公式五三：

公式五四：

在步骤四中，判定当前是否满足以及若是(即和同时满足时)，执行步骤六；否则(即和中任一个未满足，或二者都不满足时)，执行步骤五。

其中，ε₁和ε₂为0至1之间的实数，例如，ε₁＝ε₂＝0.1。

在步骤五中，根据几何变换函数F的当前估计值对当前Z_M(x_i,y_j,z_k)进行变换，将Z_M(x_i,y_j,z_k)变换后的值记为并用更新Z_M(x_i,y_j,z_k)后返回执行步骤二。

步骤六、将参数Θ的当前估计值确定为计算结果，以完成图像配准。

其中，在定义几何变换函数F(x,y,z)时，(x,y,z)是f_M(x,y,z)的像素点，而不是f_R(x,y,z)中的；同样，这个讨论同样是基于||F(x,y,z)-(x,y,z)||小于预设实数，以使f_M(F₁(x,y,z),F₂(x,y,z),F₃(x,y,z))的泰勒展开式有效；在实际中，从f_R(x,y,z)中获取像素点(x,y,z)要更合理。

优选实施例

该方法适用于一个图像到另一图像间的几何刚性变换的情况，即同一幅图中的任意两个像素间的欧氏距离变换后不变。刚体变换在实践中是最为常见的。本优选实施例采用泰勒展开式和最小二乘法，通过这样的方法，可以很好地得出参数Θ＝(α,β,γ,Δx,Δy,Δz)^T的估计值。这使得可以在3D图像采集数据大幅减少的情况下进行有效计算。

在数学上，3D图像配准问题描述如下：参考3D图像设为f_R(x,y,z)(即上文所述的参考图像)，移动3D图像设为f_M(x,y,z)(即上文所述的移动图像)。3D图像配准的主要目的是为了找到一个令f_M(F(x,y,z))尽可能逼近f_R(x,y,z)的几何变换F(x,y,z)＝(F₁(x,y,z),F₂(x,y,z),F₃(x,y,z))。对刚性变换来说，几何变换核心参数包括3D旋转、3D平移以及尺度缩放。在该优选实施例中假设没有尺度缩放，即尺度缩放系数为1，这也与大部分实际应用相符合。

假设α、β和γ分别表示沿着x轴、y轴和z轴的旋转角。(Δx,Δy,Δz)是沿着三个轴(即x轴、y轴和z轴)的平移参数，则3D刚体变换可用下面的表达式(1)来描述：

实际上，真实的图像f_R(x,y,z)和f_M(x,y,z)通常是不存在的。实际采集的图像含有噪声，由此可以建立模型如下:

在式(2)中，i＝1,2,…,n，j＝1,2,…,n，k＝1,2,…,n，其中n为正整数。例如，n可以为1000或10000，或者n也可以为无穷大整数。

其中的(x_i,y_j,z_k)∈Ω表示的是对应位置的像素；Z_R(x_i,y_j,z_k)表示实际的参考图像，而Z_M(x_i,y_j,z_k)表示实际的移动图像；ε_R(x_i,y_j,z_k)表示Z_R(x_i,y_j,z_k)中所含的噪声，而ε_M(x_i,y_j,z_k)表示Z_M(x_i,y_j,z_k)中所含的噪声。

因此，为了解决图像的配准问题，需从实际移动图像Z_M(x_i,y_j,z_k)和实际参考图像Z_R(x_i,y_j,z_k)中估计出式(1)中的(α,β,γ)和(Δx,Δy,Δz)这六个参数值，以确保Z_M(F(x_i,y_j,z_k))尽可能地逼近Z_R(x_i,y_j,z_k)，可以按照下面描述的过程来估计(α,β,γ)和(Δx,Δy,Δz)。

首先，修改几何变换F(x,y,z)＝(F₁(x,y,z),F₂(x,y,z),F₃(x,y,z))为下式表达形式：

其中，b(x,y,z)、c(x,y,z)和d(x,y,z)分别表示在x轴、y轴和z轴上的变化量，b(x,y,z)＝F₁(x,y,z)-x，c(x,y,z)＝F₂(x,y,z)-y，d(x,y,z)＝F₃(x,y,z)-z。因此，对F(x,y,z)估计，相当于对(b(x,y,z),c(x,y,z),d(x,y,z))估计。假若(b(x,y,z),c(x,y,z),d(x,y,z))很小，而且f_M(x,y,z)有在(x,y,z)上的一阶偏导数。由泰勒展开式，可以得到:

上式中，f′_Mx(x,y,z)、f′_My(x,y,z)和f′_Mz(x,y,z)分别表示f_M(x,y,z)对x、y和z的一阶偏导数，||·||表示欧几里得范数。

为了计算于(b(x,y,z),c(x,y,z),d(x,y,z))的估计需令下式的计算误差尽可能小：

f_R(x,y,z)-[f_M(x,y,z)+f′_Mx(x,y,z)b(x,y,z)+f′_My(x,y,z)c(x,y,z)+f_M'_z(x,y,z)d(x,y,z)]。

然而，事实上，f_R(x,y,z)、f_M(x,y,z)、f_M'_x(x,y,z)、f_M'_y(x,y,z)和f_M'_z(x,y,z)都是不存在的。只有式(2)中定义的有噪声的3D图像{Z_R(x_i,y_j,z_k)}和{Z_M(x_i,y_j,z_k)}是可以得到的。

在本优选实施例中，使用统计回归分析中的最小二乘法(LS)对上述变换中的六个参数进行估计，具体如下：首先，假设f′_Mx(x,y,z)、f′_My(x,y,z)和f′_Mz(x,y,z)事先已经由对应的估计函数和估计出来，具体估计方法见公式(6)。于是，b(x,y,z),c(x,y,z),d(x,y,z)中的参数Θ＝(α,β,γ,Δx,Δy,Δz)^T便可由最小化问题的解决方案来估计近似：

其中，K_h是一个以单位圆支撑的核函数，h＞0是一个带宽参数。假设参数α、β和γ都很小，那么sinα≈α，sinβ≈β，sinγ≈γ，cosα≈1，cosβ≈1，cosγ≈1。在这种条件下，由式(4)得出的Θ的最小二乘法估计值由下式计算：

Θ＝(X^TX)^-1X^TY (5)

其中X是一个n³×6设计矩阵，它的第[(i-1)n²+(j-1)n+k]行元素是由下式确定的：

此外，Y是一个n³维向量，它的第[(i-1)n²+(j-1)n+k]元素由Z_R(x_i,y_j,z_k)-Z_M(x_i,y_j,z_k)确定。

为了计算式(4)和式(5)所定义的参数估计值，可根据式(6)得到和的估计值：

核函数例如采用K_h(x,y,z)＝(1-x²)(1-y²)(1-z²)。

只有在Θ＝(α,β,γ,Δx,Δy,Δz)^T中的所有参数值很小时，几何变换F(x,y,z)才能被正确的估计。对于3D移动图像比参考图像变换较大的时候，可以使用迭代运算。由此可以通过上文结合图1描述的步骤一至步骤六完成迭代计算，获得几何变换参数，进而完成图像配准。

在该3D图像快速刚性配准方法中，通过使用3D图像灰度值建立非参数变换统计模型(例如)，突破已有方法必须满足特定的条件、需要预先设定参数等局限性。

尽管根据有限数量的实施例描述了本发明，但是受益于上面的描述，本技术领域内的技术人员明白，在由此描述的本发明的范围内，可以设想其它实施例。此外，应当注意，本说明书中使用的语言主要是为了可读性和教导的目的而选择的，而不是为了解释或者限定本发明的主题而选择的。因此，在不偏离所附权利要求书的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。对于本发明的范围，对本发明所做的公开是说明性的，而非限制性的，本发明的范围由所附权利要求书限定。

Claims

1.非特征的3D图像快速刚性配准方法，其特征在于，所述非特征的3D图像快速刚性配准方法包括：

步骤一、获得实际参考图像Z_R(x_i,y_j,z_k)和实际移动图像Z_M(x_i,y_j,z_k)；执行步骤二；

步骤二、估计下式中的和

<mrow> <msubsup> <mover> <msub> <mi>f</mi> <mi>M</mi> </msub> <mo>^</mo> </mover> <mi>x</mi> <mo>&prime;</mo> </msubsup> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>x</mi> </mrow> <mo>)</mo> </mrow> <msub> <mi>Z</mi> <mi>M</mi> </msub> <mrow> <mo>(</mo> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> </mrow> <mo>)</mo> </mrow> <msub> <mi>K</mi> <mi>h</mi> </msub> <mrow> <mo>(</mo> <mrow> <mfrac> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>x</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mi>z</mi> </mrow> <mi>h</mi> </mfrac> </mrow> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <msup> <mrow> <mo>(</mo> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>x</mi> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <msub> <mi>K</mi> <mi>h</mi> </msub> <mrow> <mo>(</mo> <mrow> <mfrac> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>x</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mi>z</mi> </mrow> <mi>h</mi> </mfrac> </mrow> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>

<mrow> <msubsup> <mover> <msub> <mi>f</mi> <mi>M</mi> </msub> <mo>^</mo> </mover> <mi>y</mi> <mo>&prime;</mo> </msubsup> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mrow> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> </mrow> <mo>)</mo> </mrow> <msub> <mi>Z</mi> <mi>M</mi> </msub> <mrow> <mo>(</mo> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> </mrow> <mo>)</mo> </mrow> <msub> <mi>K</mi> <mi>h</mi> </msub> <mrow> <mo>(</mo> <mrow> <mfrac> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>x</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mi>z</mi> </mrow> <mi>h</mi> </mfrac> </mrow> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <msup> <mrow> <mo>(</mo> <mrow> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <msub> <mi>K</mi> <mi>h</mi> </msub> <mrow> <mo>(</mo> <mrow> <mfrac> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>x</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mi>z</mi> </mrow> <mi>h</mi> </mfrac> </mrow> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>

<mrow> <msubsup> <mover> <msub> <mi>f</mi> <mi>M</mi> </msub> <mo>^</mo> </mover> <mi>z</mi> <mo>&prime;</mo> </msubsup> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mi>z</mi> </mrow> <mo>)</mo> </mrow> <msub> <mi>Z</mi> <mi>M</mi> </msub> <mrow> <mo>(</mo> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> </mrow> <mo>)</mo> </mrow> <msub> <mi>K</mi> <mi>h</mi> </msub> <mrow> <mo>(</mo> <mrow> <mfrac> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>x</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mi>z</mi> </mrow> <mi>h</mi> </mfrac> </mrow> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <msup> <mrow> <mo>(</mo> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mi>z</mi> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <msub> <mi>K</mi> <mi>h</mi> </msub> <mrow> <mo>(</mo> <mrow> <mfrac> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>x</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> </mrow> <mi>h</mi> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mi>z</mi> </mrow> <mi>h</mi> </mfrac> </mrow> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>;</mo> </mrow>

其中，和分别表示移动图像f_M(x,y,z)在(x,y,z)上的一阶偏导数f_M′_x(x,y,z)、f_M′_y(x,y,z)和f_M'_z(x,y,z)的估计函数；(x_i,y_j,z_k)∈Ω表示参考图像与移动图像对应位置的像素点，i,j,k指的是各像素点的位置下标，Z_M(x_i,y_j,z_k)表示的是实际移动图像，K_h是以单位圆支撑的核函数，K_h(x,y,z)＝(1-x²)(1-y²)(1-z²)，h＞0是一个带宽参数；执行步骤三；

步骤三、由Θ＝(X^TX)^-1X^TY获取参数Θ的当前估计值其中，Θ＝(α,β,γ,Δx,Δy,Δz)^T，

其中，X是n³×6设计矩阵，它的第[(i-1)n²+(j-1)n+k]行元素是由确定的；Y是n³维向量，它的第[(i-1)n²+(j-1)n+k]元素由Z_R(x_i,y_j,z_k)-Z_M(x_i,y_j,z_k)确定；

根据参数Θ的当前估计值确定几何变换函数F的当前估计值其中，执行步骤四；

步骤四、判定当前是否满足以及若是，执行步骤六；否则，执行步骤五；其中，ε₁和ε₂为0至1之间的实数；

步骤五、根据几何变换函数F的当前估计值对当前Z_M(x_i,y_j,z_k)进行变换，将Z_M(x_i,y_j,z_k)变换后的值记为并用更新Z_M(x_i,y_j,z_k)后返回执行步骤二；

2.根据权利要求1所述的非特征的3D图像快速刚性配准方法，其特征在于，ε₁＝ε₂＝0.1。

3.根据权利要求1所述的非特征的3D图像快速刚性配准方法，其特征在于，在定义几何变换函数F(x,y,z)时，(x,y,z)是f_M(x,y,z)的像素点；||F(x,y,z)-(x,y,z)||小于预设实数，以使f_M(F₁(x,y,z),F₂(x,y,z),F₃(x,y,z))的泰勒展开式有效。