CN103745467A - 一种基于数字体积相关法的三维图像配准方法 - Google Patents

Abstract

本发明公开了一种基于数字体积相关法的三维图像配准方法，包括如下步骤：步骤1，对样品进行第一次三维成像，得到三维图像数据G₁（x,y,z）；步骤2，刚性运动后，对样品进行第二次三维成像，获得三维图像数据G₂（x,y,z）；步骤3，对G₁（x,y,z）和G₂（x,y,z）进行数字体积相关计算，获得G₂（x,y,z）三维图像数据中部分位置点相对G₁（x,y,z）的位移值；步骤4，将步骤3得到的相对位移值代入刚性运动位移公式（1）中，采用最小二乘法计算，得到待测样品的6个刚性运动自由度；步骤5，将步骤4得到的6个刚性运动自由度代入公式（2）计算G₂（x,y,z）中各个点刚性运动后的新位置；步骤6，利用插值算法计算得到刚性运动后的新三维图像数据G₃（x,y,z）。本发明方法避免了难以控制的图像分割过程，得到的配准图像配准效果好，配准精度高。

Description

一种基于数字体积相关法的三维图像配准方法

技术领域

本发明涉及一种基于数字体积相关法的三维图像配准方法，属于图像分析技术领域。

背景技术

图像配准是图像处理的重要任务之一，主要用于将不同时间、不同视角、不同拍摄条件、甚至不同传感器获取的关于同一目标的两幅图像进行匹配的过程。图像配准的基本问题就是提出一种图像变换方法，用以校正图像的坐标和形变。把图像配准所要解决的问题抽象出共同的特性来：即在变换空间中寻找一种特定的最优的变换，达到使两幅图像在某种意义上匹配为目的。对于不同图像数据和不同应用背景的配准问题，就需要具体问题具体分析。现有的配准方法主要有基于特征的配准方法和基于区域灰度的配准方法。基于区域灰度的方法中，通常利用区域相关的方法提取控制结构，对图像的灰度依赖比较强，对图像的缩放、旋转很敏感；基于特征的方法不是直接对图像的灰度进行操作，更多的是在特征空间里提取控制结构，实现图像配准，然而特征点的选取需要用到图像分割，图像分割难以自动化，人为因素会造成较大的影响。

在很多情况下物体的变化可认为是刚性运动过程，即物体内部各位置点的相对位置未发生变化的运动过程，本专利针对的配准也限于刚性运动的配准，即刚性配准。

数字体积相关法(英文全称：Digital Volume Correlation，缩写：DVC)是由B.K.Bay等在1999年(Bay B K,Smith T S,Fyherie D P,et al.Digital volume correlation:Three-dimensional strain mapping using X-ray tomography[J].ExperimentalMechanics,1999,39(3):217～226)在传统的数字图像相关的基础上首次提出的。对变形前后采集到的物体内部的三维数字图像进行相关计算，可获得样品各位置点的位移值。利用DVC可以直接计算出样品中不同位置点的位移值u,v,w，从而可相应的计算出样品刚性运动的6个自由度，本专利据此开发了一种与传统配准方法完全不同的三维图像配准方法。

发明内容

发明目的：本发明所要解决的技术问题是提供一种基于数字体积相关法的三维图像配准方法。

发明内容：为解决上述技术问题，本发明所采用的技术方案为：

一种基于数字体积相关法的三维图像配准方法，包括如下步骤：

步骤1，对待测样品进行第一次三维成像，获得第一次三维图像数据G₁（x,y,z），将第一次三维图像作为参考图像；

步骤2，将待测样品进行刚性运动后进行第二次三维成像，获得第二次三维图像数据G₂（x,y,z），将第二次三维图像作为待配准图像；

步骤3，对G₁（x,y,z）和G₂（x,y,z）进行数字体积相关计算，获得G₂（x,y,z）三维图像数据中部分位置点相对G₁（x,y,z）的位移值；

步骤4，将步骤3得到的相对位移值代入刚性运动位移公式中，采用最小二乘法计算，得到待测样品的6个刚性运动自由度，其中，刚性运动位移公式（1）为：

$\left[\begin{array}{c}u\\ v\\ w\end{array}\right]=\left[\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right]+\left[\begin{array}{ccc}0& -w_z& w_y\\ w_z& 0& -w_x\\ -w_y& w_x& 0\end{array}\right]\·\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right];$

步骤5，将步骤4得到的6个刚性运动自由度代入公式II计算得到G₂（x,y,z）三维图像数据中各点刚性运动后的新位置，其中，公式II为：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]+\left[\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right]+\left[\begin{array}{ccc}0& {-w}_z& w_y\\ w_z& 0& {-w}_x\\ {-w}_y& w_x& 0\end{array}\right]\·\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right];$

步骤6，将步骤5得到的各点刚性运动后的新位置数据利用插值算法即可得到刚性运动后的三维图像数据G₃（x,y,z），即配准后的三维数字图像数据G₃（x,y,z）。

其中，步骤1或步骤2中，所述三维成像采用的是X射线断层照相、核磁共振断层照相、伽马射线断层照相或共聚焦显微镜中的任意一种成像技术。

其中，所述X射线断层照相包括医用X射线断层照相、工业X射线断层照相、显微X射线断层照相或纳米X射线断层照相。

其中，步骤4中，刚性运动位移公式中，u,v,w为通过数字体积相关算法得到的部分位置点或全部位置点分别在x,y,z三个方向的位移值；Δx,Δy,Δz为上述部分位置点或全部位置点离待测样品特定转动中心在x,y,z三个方向的距离；u₀,v₀,w₀,w_x,w_y,w_z是待测样品的6个刚性运动自由度。

其中，步骤5中，公式II中，（x',y',z'）是新位置点的坐标，（x,y,z）是G₂（x,y,z）中任意点的坐标。

其中，步骤3中，所述数字体积相关算法包括现存各种数字体积相关算法。

其中，步骤3中，所述部分位置点包括全部位置点和部分位置点。

其中，步骤6中，所述插值算法包括各种现存的插值算法。

有益效果：相比于现有技术，本发明的基于数字体积相关法的三维图像配准方法具有如下特点：首先，本发明方法直接利用数字体积相关法得到相对位移值，再利用最小二乘算法得到6个样品的刚性运动自由度，物理意义明确；其次，相对于基于图像特征的配准方法，本发明方法避免了难以控制的图像分割过程，得到的配准图像配准效果好，配准精度高；最后，本发明方法能够随着数字体积相关算法的发展，使其配准精度能够相应提高。

附图说明

图1为本发明基于数字体积相关法的三维图像配准方法的流程图；

图2为本发明实施例2中待测样品的X射线断层照相法获得三维图像数据G₁（x,y,z）；

图3为本发明实施例2中待测样品的X射线断层照相法获得三维图像数据G₂（x,y,z）；

图4为本发明实施例2中利用插值算法得到的刚性运动后待测样品的三维图像数据G₃（x,y,z）；

图5为本发明方法实施例2的配准效果展示图，其中，图5（a）是将G₃（x,y,z）与G₁（x,y,z）的某个截面进行图像减法运算后得到的二维截面图；图5（b）是图5（a）横线上的一维灰度图。

具体实施方式

下面结合附图和实施例，对本发明的技术方案进行详细说明，实施例用泡沫铝样品的X射线断层扫描实验来说明本发明。但本发明要求保护的范围不局限于实施例表示的范围，故凡依本发明专利申请范围所述的构造、特征及原理所做的等效变化或修饰，均包括于本发明专利申请范围内。

数字体积相关法的英文名称缩写：DVC。

刚性运动，即物体内部各位置点的相对位置未发生变化的运动过程，本发明方法针对的配准也限于刚性运动的配准，即刚性配准。

刚性配准可以用三个平动自由度u₀,v₀,w₀和三个转动自由度w_x,w_y,w_z进行描述，刚性运动所引起的物体内部每一个位置的位移值u,v,w变化满足如下的刚性运动公式：

$\left[\begin{array}{c}u\\ v\\ w\end{array}\right]=\left[\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right]+\left[\begin{array}{ccc}0& -w_z& w_y\\ w_z& 0& -w_x\\ -w_y& w_x& 0\end{array}\right]\·\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right],$ 公式中，Δx,Δy,Δz是位置点（某一点）离待测样品特定转动中心在x,y,z三个方向的距离。

实施例1：泡沫铝样品的X射线断层扫描实验，采用没有人为移动的数据对本发明配准方法进行评价：

步骤1，将泡沫铝割切成40 mm×40 mm×10 mm小块，放在工业X射线断层照相设备的旋转台上对样品进行第一次X射线断层照相成像，获得第一次三维图像数据G₁（x,y,z），将第一次三维图像数据G₁（x,y,z）作为参考图像；

步骤2，保持样品不动，连续进行第二次X射线断层照相成像，获得第二次三维图像数据G₂（x,y,z），将第二次三维图像G₂（x,y,z）作为待配准图像；

步骤3，对G₁（x,y,z）和G₂（x,y,z）进行数字体积相关计算（该数字体积相关计算参见文献Pan B, Wu D, Wang Z (2012) Internal displacement and strainmeasurement using digital volume correlation: a least-squares framework [J].Measurement Science and Technology, 23(4): 045002中的算法进行计算），得到G₂（x,y,z）三维图像数据中部分位置点（包括全部位置点和部分位置点）相对G₁（x,y,z）的位移值；

步骤4，将步骤3得到的相对位移值代入刚性运动位移公式中，采用最小二乘法计算，得到待测样品的6个刚性运动自由度，刚性运动位移公式为：

$\left[\begin{array}{c}u\\ v\\ w\end{array}\right]=\left[\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right]+\left[\begin{array}{ccc}0& -w_z& w_y\\ w_z& 0& -w_x\\ -w_y& w_x& 0\end{array}\right]\·\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right],$

其中，刚性运动位移公式中，u,v,w为通过数字体积相关算法得到的部分位置点分别在x,y,z三个方向的位移值；Δx,Δy,Δz为上述部分位置点离待测样品特定转动中心在x,y,z三个方向的距离；u₀,v₀,w₀,w_x,w_y,w_z是待测样品的6个刚性运动自由度，得到的6个刚性运动自由度为(-0.0066,0.0052,-0.1531,-1.56×10^-4,0.72×10^-4,-0.0022)。

在没有人为移动的情况下，连续两次断层成像数据的配准结果，理论上6个自由度应该均为零。但是由于在X射线断层扫描过程中样品台必不可免的平动和转动，没有办法真正回复到原来的位置，这样使样品有一定的刚性运动发生，所以其6个自由度并非全部为0。用本配准方法获得的6个自由度中有1个平动自由度和一个转动自由度较大，正是因为样品台的沿着1个方向的平动和1个方向的转动引起的；另外四个自由度较小，说明了本配准方法的高精度。

实施例2：泡沫铝样品的X射线断层扫描实验，采用有人为移动的数据对本发明配准方法进行评价：

步骤1，将泡沫铝割切成40mm×40mm×10mm小块，放在工业X射线断层照相设备的旋转台上对样品进行第一次X射线断层照相成像，获得第一次三维图像数据G₁（x,y,z），如图2所示，将第一次三维图像数据G₁（x,y,z）作为参考图像；

步骤2，将样品随意拿起来再放下，然后进行第二次X射线断层照相成像，获得第二次三维图像数据G₂（x,y,z），如图3所示，将第二次三维图像G₂（x,y,z）作为待配准图像；

步骤3，对G₁（x,y,z）和G₂（x,y,z）采用实施例1步骤3中的算法进行数字体积相关计算，得到G₂（x,y,z）三维图像数据中部分位置点（包括全部位置点和部分位置点）相对G₁（x,y,z）的位移值；

步骤4，将步骤3得到的相对位移值代入刚性运动位移公式中，采用最小二乘法计算，得到待测样品的6个刚性运动自由度，刚性运动位移公式为：

$\left[\begin{array}{c}u\\ v\\ w\end{array}\right]=\left[\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right]+\left[\begin{array}{ccc}0& -w_z& w_y\\ w_z& 0& -w_x\\ -w_y& w_x& 0\end{array}\right]\·\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right],$

其中，刚性运动位移公式（1）中，u,v,w为通过数字体积相关算法得到的部分位置点分别在x,y,z三个方向的位移值；Δx,Δy,Δz为上述部分位置点离待测样品特定转动中心在x,y,z三个方向的距离；u₀,v₀,w₀,w_x,w_y,w_z是待测样品的6个刚性运动自由度，得到的6个刚性运动自由度为(-7.64,-5.46,0.15,-0.03,0.03,0.08)；

步骤5，将步骤4得到的6个刚性运动自由度(-7.64,-5.46,0.15,-0.03,0.03,0.08)代入公式II计算三维图像数据G₂（x,y,z）中各个点刚性运动后的新位置，公式II为：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]+\left[\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right]+\left[\begin{array}{ccc}0& {-w}_z& w_y\\ w_z& 0& {-w}_x\\ {-w}_y& w_x& 0\end{array}\right]\·\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right],$

其中，公式II中，（x',y',z'）是新位置点的坐标，（x,y,z）是G₂（x,y,z）中任意点的坐标；

步骤6，将步骤5得到的各点刚性运动后的新位置利用插值算法即可得到刚性运动后的新三维图像数据G₃（x,y,z），即配准后的三维数字图像数据G₃（x,y,z），如图4所示。

将G₃（x,y,z）和G₁（x,y,z）的同一个特定截面进行图像减法运算，如图5所示，其中图5(a)为二维截面图，图5（b）是图5（a）横线上的一维灰度图，说明了本发明方法较好的配准效果。

Claims (5)

1.一种基于数字体积相关法的三维图像配准方法，其特征在于：包括如下步骤：

步骤1，对待测样品进行第一次三维成像，获得第一次三维图像数据G₁（x,y,z）；

步骤2，刚性运动后，对待测样品进行第二次三维成像，获得第二次三维图像数据G₂（x,y,z）；

步骤3，对G₁（x,y,z）和G₂（x,y,z）进行数字体积相关计算，获得G₂（x,y,z）三维图像数据中部分位置点相对G₁（x,y,z）的位移值；

步骤4，将步骤3得到的相对位移值代入刚性运动位移公式中，采用最小二乘法计算，得到待测样品的6个刚性运动自由度；

步骤5，将步骤4得到的6个刚性运动自由度代入公式II计算三维图像数据G₂（x,y,z）中各个点刚性运动后的新位置；

步骤6，将步骤5得到的各个点刚性运动后的新位置数据利用插值算法即可得到刚性运动后的三维图像数据G₃（x,y,z），即配准后的三维数字图像数据G₃（x,y,z）。

2.根据权利要求1所述的基于数字体积相关法的三维图像配准方法，其特征在于：步骤1或步骤2中，所述三维成像采用的是X射线断层照相、核磁共振断层照相、伽马射线断层照相或共聚焦显微镜中的任意一种成像技术。

3.根据权利要求2所述的基于数字体积相关法的三维图像配准方法，其特征在于：所述X射线断层照相包括医用X射线断层照相、工业X射线断层照相、显微X射线断层照相或纳米X射线断层照相。

4.根据权利要求1所述的基于数字体积相关法的三维图像配准方法，其特征在于：步骤4中，所述刚性运动位移公式为：

$\left[\begin{array}{c}u\\ v\\ w\end{array}\right]=\left[\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right]+\left[\begin{array}{ccc}0& -w_z& w_y\\ w_z& 0& -w_x\\ -w_y& w_x& 0\end{array}\right]\·\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right],$ 其中，公式中，u,v,w为通过数字体积相关算法得到的部分位置点分别在x,y,z三个方向的位移值；Δx,Δy,Δz为上述部分位置点离待测样品特定转动中心在x,y,z三个方向的距离；u₀,v₀,w₀,w_x,w_y,w_z是待测样品的6个刚性运动自由度。

5.根据权利要求1所述的基于数字体积相关法的三维图像配准方法，其特征在于：步骤5中，所述公式II为：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]+\left[\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right]+\left[\begin{array}{ccc}0& {-w}_z& w_y\\ w_z& 0& {-w}_x\\ {-w}_y& w_x& 0\end{array}\right]\·\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right],$ 其中，公式II中，（x',y',z'）是新位置点的坐标，（x,y,z）是G₂（x,y,z）中任意点的坐标。

Images