CN103310458B - 结合凸包匹配和多尺度分级策略的医学图像弹性配准方法 - Google Patents

Abstract

本发明提出了一种结合凸包匹配和多尺度分级策略的医学图像弹性配准方法，解决了大角度变换的图像配准问题，并且提高了图像配准的精度。基于凸包匹配实现点云预配准，即通过等值面提取算法提取体数据表面结构，使用点云凸包表面匹配的方法获取最优刚性配准结果；通过对图像进行多尺度分块实现弹性配准，利用多尺度迭代过程对图像做不同尺度的模糊，实现由粗至细的匹配过程；以互信息量作为相似性测度，结合三线性部分体积分布插值算法，通过迭代优化得到最优弹性变换参数。

Description

结合凸包匹配和多尺度分级策略的医学图像弹性配准方法

技术领域

本方法是一种结合凸包匹配和多尺度分级策略的医学图像弹性配准方法，可广泛应用于临床基于多模态影像的诊断治疗、病情监测和外科手术等。

背景技术

现今的医学图像获取方式包括X射线造影术、计算机断层扫描(ComputedTomography，CT)、核磁共振成像(MagneticResonanceImaging，MRI)、正电子发射断层扫描(PositronEmissionTomography，PET)等，这些成像方式可以分为解剖结构成像和功能成像两大类。解剖结构成像的特点是分辨率高，能够获得人体内脏、骨骼的形态信息；而功能成像能提供人体内脏器官的功能代谢信息，但成像分辨率低。因此，不同模态图像存在信息差异，通过采用图像融合的技术，可以将不同模态获得的图像融合，有助于提高临床诊断治疗的水平。

图像融合的关键步骤之一是将不同模态的图像进行配准。Brown早于1992年发表了较全面的关于配准技术的综述，根据变换类型的不同可以分为刚性配准和非刚性配准。刚性配准将物体视为刚体，即认为物体内部任意两点间的距离在变换前后保持不变，在医学中对应密度较高的组织结构、不易发生形变的器官。因此若变换只包括旋转、平移和缩放，则称为刚性变换。而非刚性变换则包括仿射变换、投影变换和曲线变换等多种变换方式。针对医学图像，通过不同成像方式获取的人体组织图像之间常常包含如呼吸运动和心率运动引起的不规则变化的结构，因此非刚性配准在医学图像配准中应用最为广泛。非刚性配准方法主要分为两类：基于空间变换的方法，主要有多项式函数法、薄板样条法和基函数法；基于伪物理模型方法，主要有基于偏微分方程的方法和基于B样条的形变模型。Pluim在综述中表明，互信息量作为配准的判据展现了广泛的适用性和强健的鲁棒性，即一幅图像中的灰度区域在另一幅图像中对应的区域也应具有相近的灰度分布。互信息量是基于信息论，对图像灰度进行概率统计的一种相似度测量法，当两幅图像完全配准时，它们的互信息量达到最大值。

由于医学图像变换较为复杂，目前图像融合采用的配准方法主要是弹性配准方法或刚性与弹性配准相结合的方法。这些传统方法的缺陷主要包括：1)当图像存在大角度偏转或者大尺度位移时，无论是基于空间变换还是基于伪物理模型方法的非刚性变换方法都会有较大的误差，甚至在特殊情况下难以实现配准。2)弹性配准的精度亦是随着图像的质量变化，是难以准确控制的因素。在配准过程中，过高的精度会使得图像配准过程耗时长，效率低，而过低的精度会使得配准的结果误差大。

发明内容

本发明提出了一种结合凸包匹配和多尺度分级策略的医学图像弹性配准方法，解决了大角度变换的图像配准问题，并且提高了图像配准的精度。

本发明的结合凸包匹配和多尺度分级策略的医学图像弹性配准方法，包括以下步骤：

(1)使用等值面提取算法分割出体数据表面结构，通过均匀采样的方式得到图像的点云数据，包括原始点云和目标点云；

(2)对所述的原始点云和目标点云分别计算其三维凸包结构及其凸包表面的有向三角形集合；

(3)对原始点云的任意一个有向三角形，与目标点云中的每一潜在有向三角形进行配对，计算出两个三角形之间的坐标变换参数，作为刚性配准的候选参数；

(4)以变换后的原始点云与目标点云间相似程度最大为目标，使用优化算法筛选获得最优刚性变换参数；

(5)基于获得的最优刚性预配准的结果，对源图像和参考图像做高斯模糊，设定一个模糊核序列，使得该序列中模糊核尺度连续降低，将此序列作为外部迭代方向；

(6)将源图像和参考图像以逐级分块的方式分别划分为若干个子块，并将该两幅图像中具有相同空间位置关系的子块作为对应子块，划分子块数随迭代层数逐层上升；

(7)计算该两幅图像中对应子块间的刚性变换参数，即通过迭代坐标变换矩阵参数，使得对应子块间的相似性测度达到最优值；

(8)使用各组对应子块的刚性变换矩阵对源图像子块进行形变，子块各自的局部形变组成的全局形变对于整幅图像表现为弹性形变；

(9)计算形变之后的整幅图像与参考图像间的相似性测度，以互信息量作为相似性测度；

(10)若此次迭代相似性测度大于上次迭代结果，则重复步骤(7)-(9)，直至相似性测度不再增加；

(11)降低高斯模糊核的尺度，重复步骤(5)-(10)，直至模糊核尺度降到最低，最终得到各个子块的刚性变换矩阵，从每个子块选取控制点，用薄板样条函数表示图像的全局变换。

本方法的优点主要体现在：1)适用性广，可用于存在复杂模型变换的图像，包括大尺度刚性变换、弹性变换或者是两者的组合；2)效率高。预配准模块中使用点云凸包配准方法，基于凸包表面矢量三角形的特性，减少了匹配计算次数，计算复杂度低；3)配准结果精度和鲁棒性高。采用多尺度的方法，由粗至细的对图像进行多层次的配准。

附图说明

图1多尺度弹性配准流程图。

具体实施方式

下面结合附图对本发明做进一步介绍。

如图1所示，基于凸包匹配实现点云预配准，即通过等值面提取算法提取体数据表面结构，使用点云凸包表面匹配的方法获取最优刚性配准结果；通过对图像进行多尺度分块实现弹性配准，利用多尺度迭代过程对图像做不同尺度的模糊，实现由粗至细的匹配过程；以互信息量作为相似性测度，结合三线性部分体积分布插值算法，通过迭代优化得到最优弹性变换参数。

本方法目的是对两幅图像进行配准，选择图像I_x作为待配准的原始图像，I_y作为参考的目标图像。本专利通过三个部分描述算法：体数据表面提取、刚性预配准和多尺度弹性配准，这三部分在流程上是连续的，详细介绍如下。

1、基于等值面提取算法的体数据表面提取

如果把体数据看成是三维空间内关于某种物理属性的采样集合，非采样点上的值可用邻近采样点插值来估计，则将该空间区域内具有相同值的点的集合定义为一个或多个曲面，称之为等值面。由于不同的组织具有不同的密度，可选定适当的值来定义等值面，该等值面表示不同物质的交界。

等值面是空间中所有具有相同密度范围的点的集合，它可以表示成{(x,y,z),f(x,y,z)＝c}，这里c是常数，通过人为设定该值选取出体数据的表面。MC(MarchingCube)算法将八个相邻的采样点包含的立方体区域定义为一个集合，当体素的部分角点值大于c，而另一部分角点小于c，则表明该体素与等值面相交。当把所有与等值面相交的体素都找到，并求出相应的相交面后，等值面也就求出来了。

这里通过定义适当的c值，可以分别提取出原始图像和目标图像的体数据轮廓点，即两组点云，原始点云和目标点云

2、基于凸包匹配的刚性预配准

这部分的目的是对图像进行全局刚性配准，忽略图像内部的变换，最终得到初始刚性变换矩阵。如图1所示，详细过程如下：

1)对于原始点云P_X和目标点云P_Y，计算其三维凸包结构及其有向三角形集合F(P_X)和F(P_Y)。

2)对原始点云的任意一个有向三角形F(P_X)_i，与目标点云的有向三角形F(P_Y)_j进行配对，估计出两个三角形之间的坐标变换参数其中R为旋转矩阵、为平移向量、s为尺度缩放因子。

3)使F(P_X)_i经过变换后与F(P_Y)_j之间的相似度最大，即对应点云之间的欧氏距离最小，该过程采用Levenberg-Marquarat(LM)算法，目标函数可以概括为：

$T{(R,\stackrel{\→}{t},s)}_{i,j}^{*}=\arg \min d(T\·F{\left(P_X\right)}_i,F{\left(P_Y\right)}_j)$

其中表示F(P_X)_i与F(P_Y)_j配对后对应的最优变换参数。将作用于原始点云，可得变换后的点云与目标点云间的欧氏距离。

4)任意三角形F(P_X)_i与目标三角形集合之间存在最优一致关系，使得对应点云间的欧氏距离在所有配对中最小。

此时F(P_X)_i与其最优配对间的变换矩阵即为F(P_Y)_j对应的最优变换对每一个三角形F(P_Y)_j与F(P_Y)进行配对，在每组配对结果中取最优的一致性关系，最终获得变换矩阵如下：

$T{(R,\stackrel{\→}{t},s)}_i^{*}=\arg \max I\left(T\right(R,\stackrel{\→}{t},s)\·P_X,P_Y)$

在刚性配准结束后，将得到的变换矩阵作用于待配准的原始图像I_x，刚性变换矩阵包括欧式变换以及尺度缩放两部分。欧氏变换受到旋转矩阵与平移向量共同约束，旋转矩阵R由绕三轴旋转α,β,γ3个参数组成，平移向量则有t_x,t_y,t_z3个参数，即欧氏变换由6个参数限定。尺度缩放变换由缩放因子s限定。因此三维坐标系下的刚性变换矩阵为7个参数约束下的坐标变换矩阵，即：

$T(R,\stackrel{\→}{t},s)=\left[\begin{array}{cccc}\cos \mathrm{\β}\cos \mathrm{\γ}\·s& \sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}& \cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}+\sin \mathrm{\α}\sin \mathrm{\γ}& t_x\\ \cos \mathrm{\β}\sin \mathrm{\γ}& (\sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}+\cos \mathrm{\α}\cos \mathrm{\γ})\·s& \cos \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}-\sin \mathrm{\α}\cos \mathrm{\γ}& t_y\\ -\sin \mathrm{\β}& \sin \mathrm{\α}\cos \mathrm{\β}& \cos \mathrm{\α}\cos \mathrm{\β}\·s& t_z\\ 0& 0& 0& 1\·s\end{array}\right]$

任意点的坐标变换可表示为：

3、多尺度弹性配准过程

多尺度弹性配准迭代过程可分为外部迭代和内部迭代。其中外部迭代为实现多尺度模糊和分块的过程，内部迭代描述了通过优化变换矩阵参数寻找最优配准的过程。详细描述如下：

1)多尺度模糊迭代优化

多尺度体现在两方面：a)医学图像中存在着不同尺度的人体结构，如血管、肌肉等。这些结构由于其物理特性具有不同的弹性变换特性，对图像进行模糊可以忽略图像中的细节。因此，若对图像做一系列高斯模糊，并且逐次降低模糊核，可以实现先对大尺度的结构进行配准，再对小尺度的结构进行配准，依次得到更精确的配准结果。b)在配准过程中对图像进行多尺度分块，实现局部配准，分块的精度随着模糊精度的变化而改变。多尺度分块采用了分治法的概念，通过将规模大的图像形变问题分解为若干个规模较小模式相同的子问题，合并子问题的解可以得到最终整幅图像的弹性变换矩阵。在弹性配准中降低了对形变模型强健性的要求。

本方法采用高斯模糊法对图像进行模糊处理，高斯模糊作用于像素后，分布不为零的像素组成的卷积矩阵与原始图像进行卷积变换。每个像素的值都是周围相邻像素值的加权平均。原始像素的值有最大的高斯分布值，因而有最大的权重，相邻像素随着距离原始像素越来越远，其权重也越来越小，这样进行模糊处理比其它的均衡模糊滤波器更高地保留了边缘信息。设定好每层迭代的模糊尺度，即模糊半径r，如{G₁,G₂,G₃}(G₁＞G₂＞G₃)，而后按照设定的顺序进行迭代。

前面提到的点云凸包匹配是对全局的刚性配准，而医学图像中人体组织常常存在局部的弹性变换，因此无法使用一个变换矩阵来描述整幅图像的变换，因此需要对图像的形变进行局部描述。本方法是对图像进行逐级分块，不同的区域进行不同的形变描述，而分块的数量决定了弹性变换的精度。理论上分块越多弹性形变精度越高，但是由于分块减少了图像的信息，例如对于基于灰度统计的互信息量，越小的分块可能会导致统计概率的波动，从而出现局部极值，针对此问题上一步对图像做模糊处理可以降低这种波动。

多尺度分块策略是对源图像和参考图像以逐级分块的方式分别划分为若干个子块，并将两幅图像中具有相同空间位置关系的子块作为对应子块。迭代过程中由粗至细、逐层的增加分块数量，直至迭代结束，得到精度最高配准结果。

2)本方法采用互信息量作为图像迭代优化的相似性测度，互信息量的表达式为：

I(T(I_x),I_y)＝H(T(I_x))+H(I_y)-H(T(I_x),I_y)

其中H(T(I_x))和H(I_y)分别是初始变换后的原始图像T(I_x)和参考图像I_y的香农熵(Shannonentropy)表示，它的值表面一幅图像中所含信息量的大小，其表达式为p_i为像素i出现的概率。H(T(I_x),I_y)为两幅图像的联合熵，可通过两幅图像的联合直方图求出联合概率分布p(i,j)，则可由求出联合熵。

3)插值策略和优化算法。待配准的源图像做了变换之后，像素点产生了位移，可能坐标位置不在整数像素上，因此需要用灰度插值的方法对像素进行估计。最近领域插值算法由于该点的值只由距该点最近的一点决定，不能达到亚像素级的精度。三线性插值算法可能会在参考图中引入新的灰度值，使得当图像发生微小变换时，参考图的灰度概率分布产生很大的变化。本方法使用三线性部分体积分布插值算法对变换后的图像进行插值，浮动图中的任意一点的灰度对联合直方图的贡献是由参考图中的点周围最近邻的八个点取与三线性插值算法相同的权重加权而得。因此变换矩阵的参数发生变化时，联合直方图的变化平缓。

优化算法是寻找最优变换参数，得到最大相似性测度。图像配准在本质上是一个多参数优化的问题，即寻找互相关系数最大、联合熵最小或互信息量最大时的几个空间变换参数值。这是一个寻找多维变量最大值的问题，常用的有Powell算法，这是一种直接搜索法，将整个搜索的过程分为几个阶段。该算法先给出n个正交的初始搜索方向第一轮寻优循环先沿这一系列方向做出一维极值搜索，即在每一个搜索方向求得该方向的近似最佳因子t，使得新的设计逼近该方向的最佳值：

$X_k^{\left(1\right)}=X_{k-1}^{\left(1\right)}+t\·S_k$

完成第一次循环后，得到的设计变量与循环开始前的设计变量之差所指方向代替原搜索方向S₁，进行第二轮循环的搜索，循环结束后，同样得出循环开始前后的设计变量之差，代替原搜索方向S₂，进行第三轮循环，如此循环下去直至满足终止的条件。对于图像变换寻找最优参数问题，包括绕三轴旋转α,β,γ3个参数、平移向量t_x,t_y,t_z3个参数和尺度缩放因子s七个参数。按照设定的顺序(t_x,t_y,γ,t_z,α,β,s)，以最大互信息量为方向，寻找此问题的多维极值。

Claims (3)

1.结合凸包匹配和多尺度分级策略的医学图像弹性配准方法，其特征在于，包括以下步骤：

(1)分割出体数据表面结构，通过均匀采样的方式得到图像的点云数据，点云数据包括原始点云和目标点云；

(2)对所述的原始点云和目标点云分别计算其三维凸包结构及其凸包表面的有向三角形集合；

(3)对原始点云的任意一个有向三角形，与目标点云中的每一潜在有向三角形进行配对，计算出两个三角形之间的坐标变换参数，作为刚性预配准的候选参数；

(4)以变换后的原始点云与目标点云间相似程度最大为目标，使用优化算法筛选获得最优刚性变换参数,获得最优刚性预配准结果；

(5)基于获得的最优刚性预配准的结果，对源图像和参考图像做高斯模糊，设定一个模糊核序列，使得该序列中模糊核尺度连续降低，将此序列作为外部迭代方向；

(6)将源图像和参考图像以逐级分块的方式分别划分为若干个子块，并将该两幅图像中具有相同空间位置关系的子块作为对应子块，划分子块数随迭代层数逐层上升；

(7)计算该两幅图像中对应子块间的刚性变换参数，即通过迭代坐标变换矩阵参数，使得对应子块间的相似性测度达到最优值；

(8)使用各组对应子块的刚性变换矩阵对源图像子块进行形变，子块各自的局部形变组成的全局形变对于整幅图像表现为弹性形变；

(9)计算形变之后的整幅图像与参考图像间的相似性测度，以互信息量作为相似性测度；

(10)若此次迭代相似性测度大于上次迭代结果，则重复步骤(7)-(9)，直至相似性测度不再增加；

(11)降低高斯模糊核的尺度，重复步骤(5)-(10)，直至模糊核尺度降到最低，最终得到各个子块的刚性变换矩阵，从每个子块选取控制点，用薄板样条函数表示图像的全局变换。

2.如权利要求1所述的结合凸包匹配和多尺度分级策略的医学图像弹性配准方法，其特征在于，步骤(1)中采用等值面提取算法分割体数据表面结构。

3.如权利要求1所述的结合凸包匹配和多尺度分级策略的医学图像弹性配准方法，其特征在于，采用下述方法进行刚性预配准：

1)对于原始点云P_X和目标点云P_Y，计算其三维凸包结构及其有向三角形集合F(P_X)和F(P_Y)；

2)对原始点云的任意一个有向三角形F(P_X)_i，与目标点云的有向三角形F(P_Y)_j进行配对，估计出两个三角形之间的坐标变换参数其中R为旋转矩阵、为平移向量、s为尺度缩放因子；

3)使F(P_X)_i经过变换后与F(P_Y)_j之间的相似度最大，即对应点云之间的欧氏距离最小，目标函数概括为：

$T{(R,\stackrel{\→}{t},s)}_{i,j}^{*}=\arg \min d(T\·F{\left(P_X\right)}_i,F{\left(P_Y\right)}_j)$

其中表示F(P_X)_i与F(P_Y)_j配对后对应的最优变换参数，将作用于原始点云，得到变换后的点云与目标点云间的欧氏距离；

4)任意三角形F(P_X)_i与目标三角形集合之间存在最优一致关系，使得对应点云间的欧氏距离在所有配对中最小；此时F(P_X)_i与其最优配对间的变换矩阵即为F(P_Y)j对应的最优变换对每一个三角形F(P_Y)_j与F(P_Y)进行配对，在每组配对结果中取最优的一致性关系，最终获得变换矩阵如下:

$T{(R,\stackrel{\→}{t},s)}_i^{*}=\arg \max I\left(T\right(R,\stackrel{\→}{t},s)\·P_X,P_Y).$