本发明涉及图像融合和图像序列监测技术领域，具体涉及一种非特征的3d图像快速刚性配准方法。
背景技术：
3d图像在应用中越来越受欢迎。例如在mri(磁共振成像)中，人们习惯于从3d物体中获取一组3d图像来研究3d物体(如：患者的头部)的生物机制。为了合理地去比较同一个对象在不同时刻两个3d图像的差异，需要事先对其进行配准(ir)。图像配准技术广泛的应用在医疗成像，目标跟踪，遥感技术，指纹和脸部识别，图像压缩和视频增强等场合。
文献中存在的大部分的2d图像配准方法大致可以分为基于特征和基于灰度方法的两类。基于特征的方法首先要考虑选取两张图像的两个特征集，然后去寻找能够很好的匹配两个特征集的几何变换f。常用的特征包括角点、lbp、surf、质心或模板，还有图像灰度函数得到的退化像素。由于特征提取通常是一个耗时多、随意性强且具有挑战性的任务，最近ir研究更多关注直接对两幅图像灰度值的变换f，包括基于参数的转换，和更灵活的非参数的转换。
一些2d图像配准方法被引入到3dslicer软件包(http://slicer.org/)中从而推广到3d情况下。最近比较流行的icp算法是基于点云数据集，点云数据是3d扫描设备产生的，主要用来代表物体的外表面形状，而现实中很多的3d图像是存放整个图像的各个位置的灰度信息(比如mri图像)，基于点云的配准算法不适用这些情况。
综上可知，现有的一些配准方法在用于3d图像配准时，往往存在配准结果不准确的问题。
技术实现要素：
在下文中给出了关于本发明的简要概述，以便提供关于本发明的某些方面的基本理解。应当理解，这个概述并不是关于本发明的穷举性概述。它并不是意图确定本发明的关键或重要部分，也不是意图限定本发明的范围。其目的仅仅是以简化的形式给出某些概念，以此作为稍后论述的更详细描述的前序。
鉴于此，本发明提出一种非特征的3d图像快速刚性配准方法，以至少解决现有配准方法在用于3d图像配准时存在配准结果不准确的问题。
根据本发明的一个方面，提供了一种非特征的3d图像快速刚性配准方法，所述非特征的3d图像快速刚性配准方法包括：步骤一、获得实际参考图像zr(xi,yj,zk)和实际移动图像zm(xi,yj,zk)；执行步骤二；步骤二、估计下式中的
其中，
进一步地，ε1＝ε2＝0.1。
进一步地，在定义几何变换函数f(x,y,z)时，(x,y,z)是fm(x,y,z)的像素点；||f(x,y,z)-(x,y,z)||小于预设实数，以使fm(f1(x,y,z),f2(x,y,z),f3(x,y,z))的泰勒展开式有效。
现有技术中，大都基于特征点和针对特定的约束条件，带来了特征选择耗时多，随机性强，而且约束条件使用不灵活等问题。相比之下，本发明的非特征的3d图像快速刚性配准方法针对这些问题，提出直接使用图像灰度值的无特征3d刚性配准方法，该方法使用泰勒展开式和最小二乘法直接计算待配准图像的变换参数，并且使用较少的数据点完成快速的配准。本发明的算法能够获得较高的精度，并且在数据减少时仍可以有效计算；而且，该方法在需要使用数据压缩的应用里是非常实用的，提出的方法适合于应用在当数据量较大或者需要快速进行图像配准场合。
通过以下结合附图对本发明的最佳实施例的详细说明，本发明的这些以及其他优点将更加明显。
附图说明
本发明可以通过参考下文中结合附图所给出的描述而得到更好的理解，其中在所有附图中使用了相同或相似的附图标记来表示相同或者相似的部件。所述附图连同下面的详细说明一起包含在本说明书中并且形成本说明书的一部分，而且用来进一步举例说明本发明的优选实施例和解释本发明的原理和优点。在附图中：
图1是示意性地示出本发明的非特征的3d图像快速刚性配准方法的一个示例的结构图。
本领域技术人员应当理解，附图中的元件仅仅是为了简单和清楚起见而示出的，而且不一定是按比例绘制的。例如，附图中某些元件的尺寸可能相对于其他元件放大了，以便有助于提高对本发明实施例的理解。
具体实施方式
在下文中将结合附图对本发明的示范性实施例进行描述。为了清楚和简明起见，在说明书中并未描述实际实施方式的所有特征。然而，应该了解，在开发任何这种实际实施例的过程中必须做出很多特定于实施方式的决定，以便实现开发人员的具体目标，例如，符合与系统及业务相关的那些限制条件，并且这些限制条件可能会随着实施方式的不同而有所改变。此外，还应该了解，虽然开发工作有可能是非常复杂和费时的，但对得益于本公开内容的本领域技术人员来说，这种开发工作仅仅是例行的任务。
在此，还需要说明的一点是，为了避免因不必要的细节而模糊了本发明，在附图中仅仅示出了与根据本发明的方案密切相关的装置结构和/或处理步骤，而省略了与本发明关系不大的其他细节。
本发明的实施例提供了一种非特征的3d图像快速刚性配准方法，所述非特征的3d图像快速刚性配准方法包括：步骤一、获得实际参考图像zr(xi,yj,zk)和实际移动图像zm(xi,yj,zk)；执行步骤二；步骤二、估计下式中的
其中，
图1给出了本发明的非特征的3d图像快速刚性配准方法的一个示例性处理。
如图1所示，该方法开始后，首先执行步骤一。
在步骤一中，获得实际参考图像zr(xi,yj,zk)和实际移动图像zm(xi,yj,zk)，然后，执行步骤二。
在步骤二中，估计公式一一至公式一三中的
公式一一：
公式一二：
公式一三：
其中，
在步骤三中，由公式二获取参数θ的当前估计值
公式二：θ＝(xtx)-1xty。
其中，θ＝(α，β，γ，δx，δy，δz)t，
其中，在参数θ中，α、β和γ分别表示沿x轴、y轴和z轴的旋转角度，而δx、δy和δz分别表示沿着x轴、y轴和z轴这三个轴的平移量。
其中，x是n3×6设计矩阵，它的第[(i-1)n2+(j-1)n+k]行元素是由公式三确定的。
公式三：
此外，y是n3维向量，它的第[(i-1)n2+(j-1)n+k]元素由zr(xi,yj,zk)-zm(xi,yj,zk)确定
这样，可以根据参数θ的当前估计值确定几何变换函数f的当前估计值
公式四一：f(x,y,z)＝(f1(x,y,z),f2(x,y,z),f3(x,y,z))
公式四二：f1(x,y,z)＝x(cosαcosβ)+y(sinαcosβ)-z(sinβ)+δx
公式四三：
f2(x,y,z)＝x(cosαsinβsinγ-sinαcosγ)+y(sinαsinβsinγ+cosαcosγ)+z(cosβsinγ)+δy
公式四四：
f3(x,y,z)＝x(cosαsinβcosγ+sinαsinγ)+y(sinαsinβcosγ-cosαsinγ)+z(cosβcosγ)+δz
公式五一：
公式五二：
公式五三：
公式五四：
在步骤四中，判定当前是否满足
其中，ε1和ε2为0至1之间的实数，例如，ε1＝ε2＝0.1。
在步骤五中，根据几何变换函数f的当前估计值
步骤六、将参数θ的当前估计值
其中，在定义几何变换函数f(x,y,z)时，(x,y,z)是fm(x,y,z)的像素点，而不是fr(x,y,z)中的；同样，这个讨论同样是基于||f(x,y,z)-(x,y,z)||小于预设实数，以使fm(f1(x,y,z),f2(x,y,z),f3(x,y,z))的泰勒展开式有效；在实际中，从fr(x,y,z)中获取像素点(x,y,z)要更合理。
优选实施例
该方法适用于一个图像到另一图像间的几何刚性变换的情况，即同一幅图中的任意两个像素间的欧氏距离变换后不变。刚体变换在实践中是最为常见的。本优选实施例采用泰勒展开式和最小二乘法，通过这样的方法，可以很好地得出参数θ＝(α,β,γ,δx,δy,δz)t的估计值。这使得可以在3d图像采集数据大幅减少的情况下进行有效计算。
在数学上，3d图像配准问题描述如下：参考3d图像设为fr(x,y,z)(即上文所述的参考图像)，移动3d图像设为fm(x,y,z)(即上文所述的移动图像)。3d图像配准的主要目的是为了找到一个令fm(f(x,y,z))尽可能逼近fr(x,y,z)的几何变换f(x,y,z)＝(f1(x,y,z),f2(x,y,z),f3(x,y,z))。对刚性变换来说，几何变换核心参数包括3d旋转、3d平移以及尺度缩放。在该优选实施例中假设没有尺度缩放，即尺度缩放系数为1，这也与大部分实际应用相符合。
假设α、β和γ分别表示沿着x轴、y轴和z轴的旋转角。(δx,δy,δz)是沿着三个轴(即x轴、y轴和z轴)的平移参数，则3d刚体变换可用下面的表达式(1)来描述：
实际上，真实的图像fr(x,y,z)和fm(x,y,z)通常是不存在的。实际采集的图像含有噪声，由此可以建立模型如下:
在式(2)中，i＝1,2,…,n，j＝1,2,…,n，k＝1,2,…,n，其中n为正整数。例如，n可以为1000或10000，或者n也可以为无穷大整数。
其中的(xi,yj,zk)∈ω表示的是对应位置的像素；zr(xi,yj,zk)表示实际的参考图像，而zm(xi,yj,zk)表示实际的移动图像；εr(xi,yj,zk)表示zr(xi,yj,zk)中所含的噪声，而εm(xi,yj,zk)表示zm(xi,yj,zk)中所含的噪声。
因此，为了解决图像的配准问题，需从实际移动图像zm(xi,yj,zk)和实际参考图像zr(xi,yj,zk)中估计出式(1)中的(α,β,γ)和(δx,δy,δz)这六个参数值，以确保zm(f(xi,yj,zk))尽可能地逼近zr(xi,yj,zk)，可以按照下面描述的过程来估计(α,β,γ)和(δx,δy,δz)。
首先，修改几何变换f(x,y,z)＝(f1(x,y,z),f2(x,y,z),f3(x,y,z))为下式表达形式：
其中，b(x,y,z)、c(x,y,z)和d(x,y,z)分别表示在x轴、y轴和z轴上的变化量，b(x,y,z)＝f1(x,y,z)-x，c(x,y,z)＝f2(x,y,z)-y，d(x,y,z)＝f3(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))很小，而且fm(x,y,z)有在(x,y,z)上的一阶偏导数。由泰勒展开式，可以得到:
上式中，f′mx(x,y,z)、f′my(x,y,z)和f′mz(x,y,z)分别表示fm(x,y,z)对x、y和z的一阶偏导数，||·||表示欧几里得范数。
为了计算于(b(x,y,z),c(x,y,z),d(x,y,z))的估计
fr(x,y,z)-[fm(x,y,z)+f′mx(x,y,z)b(x,y,z)+f′my(x,y,z)c(x,y,z)+fm'z(x,y,z)d(x,y,z)]。
然而，事实上，fr(x,y,z)、fm(x,y,z)、fm'x(x,y,z)、fm'y(x,y,z)和fm'z(x,y,z)都是不存在的。只有式(2)中定义的有噪声的3d图像{zr(xi,yj,zk)}和{zm(xi,yj,zk)}是可以得到的。
在本优选实施例中，使用统计回归分析中的最小二乘法(ls)对上述变换中的六个参数进行估计，具体如下：首先，假设f′mx(x,y,z)、f′my(x,y,z)和f′mz(x,y,z)事先已经由对应的估计函数
其中，kh是一个以单位圆支撑的核函数，h＞0是一个带宽参数。假设参数α、β和γ都很小，那么sinα≈α，sinβ≈β，sinγ≈γ，cosα≈1，cosβ≈1，cosγ≈1。在这种条件下，由式(4)得出的θ的最小二乘法估计值由下式计算：
θ＝(xtx)-1xty(5)
其中x是一个n3×6设计矩阵，它的第[(i-1)n2+(j-1)n+k]行元素是由下式确定的：
此外，y是一个n3维向量，它的第[(i-1)n2+(j-1)n+k]元素由zr(xi,yj,zk)-zm(xi,yj,zk)确定。
为了计算式(4)和式(5)所定义的参数估计值，可根据式(6)得到
核函数例如采用kh(x,y,z)＝(1-x2)(1-y2)(1-z2)。
只有在θ＝(α,β,γ,δx,δy,δz)t中的所有参数值很小时，几何变换f(x,y,z)才能被正确的估计。对于3d移动图像比参考图像变换较大的时候，可以使用迭代运算。由此可以通过上文结合图1描述的步骤一至步骤六完成迭代计算，获得几何变换参数，进而完成图像配准。
在该3d图像快速刚性配准方法中，通过使用3d图像灰度值建立非参数变换统计模型(例如
尽管根据有限数量的实施例描述了本发明，但是受益于上面的描述，本技术领域内的技术人员明白，在由此描述的本发明的范围内，可以设想其它实施例。此外，应当注意，本说明书中使用的语言主要是为了可读性和教导的目的而选择的，而不是为了解释或者限定本发明的主题而选择的。因此，在不偏离所附权利要求书的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。对于本发明的范围，对本发明所做的公开是说明性的，而非限制性的，本发明的范围由所附权利要求书限定。