CN102169578B - 一种基于有限元模型的非刚性医学图像配准方法 - Google Patents

Abstract

本发明属于医学图像处理领域，涉及一种基于有限元模型非刚性医学图像配准方法。本发明所涉及配准方法在配准过程中采用分级配准的思想，先对浮动图像和参考图像进行全局配准，再对全局配准图像和参考图像进行有限元建模，最后计算与该模型相关的能量项，包内部能量项和驱动模型发生形变的外部能量项。在内部能量项和外部能量项的共同控制下模型的整体总能量达到最小时，模型收敛于平衡位置，实现图像局部精确配准，最终实现非刚性图像的配准。实验表明本方法能够能很好的模拟非刚性图像的局部形变，从而提高配准的精确度及速度、增强配准的鲁棒性，为医学图像融合奠定一个良好的基础。

Description

一种基于有限元模型的非刚性医学图像配准方法

技术领域

本发明涉及一种图像处理技术领域的方法，具体是一种基于有限元模型的非刚性医学图像配准方法。 

背景技术

医学图像配准是医学图像融合的基础，也是医学图像分析一项重要的技术，其主要通过寻找一种空间变换，使得两幅医学图像上的对应点达到空间位置或者解剖点上的一致性。在医学诊断过程中，由于存在不同模式图像表现不同性质的物理机制、患者的移动、成像参数的变化、不同成像设备间空间分辨率不相同等现实问题，所以单单凭借医生手动将两张相同或不同模式的图像在空间上作配准会受到很多局限，且常带有较大的主观性，其可靠性往往不高，不可避免地会产生误差。特别是在定向放射外科和心脏手术可视化等应用领域，对于图像配准的精度要求很高，使得医学图像配准成为一项必要而有相当困难的任务，诸如需要进一步解决刚性变换或者放射变换等全局配准方法通常无法应付局部变换的缺点。为此本文提出了一种基于有限元模型的非刚性医学图像配准方法。 

医学图像配准方法的步骤主要包括：确定待配准图像与目标图像的一种空间变换；确定经过空间变换后的图像与目标图像的相似性测度；寻找一种参数优化策略使两幅图像的相似性测度达到最大。现有配准方法主要分为两大类：基于图像角点医学图像配准方法与基于医学图像灰度信息的配准方法。两者技术缺陷主要为： 

基于图像特征的医学图像配准方法的技术缺陷为：它需要对图像进行分割来提取图像的特征，由于非刚性组织的结构很复杂，有些分界面不是很明显，通常需要人工预选定特征，这样会费时也费力而且配准的精度受分割精度影响，一般很难自动完成，使得配准时间过程、速度慢、配准精确不高。 

基于图像灰度信息的医学图像配准方法的技术缺陷为：它不需要对图像进行分割处理，直接对整幅图像进行运算，会造成配准的速度较慢、配准时间长、鲁棒性差。 

非刚性医学图像相对刚性医学图像来讲，由于其形变量较大较复杂，现有的配准算法不能很好的模拟其形变并克服配准时消耗时间多、精确度低、鲁棒性不强等缺点。 

发明内容

本发明的目的在于提供一种能够提高配准的精确度及速度，增强配准鲁棒性的基于有限元模型的非刚性医学图像配准方法。 

技术解决方案： 

本发明一种基于有限元模型非刚性医学图像的配准方法，方法步骤如下：1)首先通过图像的一阶矩寻找浮动图像和参考图像这两幅待配准图像的质心；再次，通过图像的二阶中心矩寻找两幅待配准图像的主轴和坐标系的夹角；最后，通过平移和旋转使得两幅待配准图像的质心和主轴对齐，得到旋转分量，实现两幅图像的全局配准。 

2)然后对全局配准图像及参考图像这两幅待配准图像进行网格划分，首先，对全局配准图像和参考图像分别进行高斯平滑滤波，消除图像成像过程中生成的噪声；然后分别对平滑处理后的两幅待配准图像利用Harris算子提取图像的角点，对其角点进行排序，其横坐标中最小值、最大值，取其纵坐标中最小值、最大值并记录下来，这样就形成四个点，连结四个角点形成一个四边形，使得四边形包含所有的角点。本发明将四边形进行划分成为横纵坐标相等的区域，将形成规则的四边形，接着按照正对角线连接四边形单元的对角线形成三角形网格，对三角形网格进行编号1、2、......、M，M为大于1的整数。 

3)建立数学模型：将待配准图像的形变过程看作弹性的伸缩过程，利用线性弹性模型中的最小势能原理约束每个匹配后的形变，外力为两幅图像的灰度差值，当外力为零时有限元单元的势能到最小，图像完成精确配准。 

4)建立有限元方程：在对步骤3)进行求解的过程中。由于对整幅图像进行建模，使得方程求解过程比较复杂，所以将浮动图像和参考图像进行离散化，将图像感兴趣的区域划分成规则形的三角形网格，将整幅图像的形变看作为三角形单元形变的总和。如步骤2)所述。将步骤2)全局配准图像和参考图像这两幅待配准图像网格单元建立有限元方程，首先，对两幅待配准图像的网格单元进行匹配，将匹配后的网格单元，进行编号1、2、......、N，N为大于1的整数：然后图像的总势能看作为三角形势能的总和，建立有限元方程： 

其中，σ^t表示弹性体的应力，ε表示弹性体的应变，∫_Ωσ^tεdΩ表示弹性系统的形变势能，Ω表示二维图像区域，∫_ΩF_iudΩ表示弹性系统的外力势能，F_i表示作用于非刚性医学图像的外力，u表示图像的形变位移，m表示三角形单元的数目，i表示第i个三角形； 

5)求解方程：求解步骤4)中的有限元方程，在每个单元内选择至少一个节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，通过插值函数计算出各个单元内场函数的近似解；利用三角形单元的灰度差值的平方和作为外力，当总势能最小时，外力为零，来求解方程，得到全局配准图像的偏移量，实现图像的局部精确配准。 

本发明将待配准的浮动图像和参考图像的质心和主轴对齐，实现两幅待配准图像的全局配准；然后根据有限分析中的变分原理将待配准图像划分成为不同的矩形区域，利用线性弹性模型最小势能原理约束图像的局部形变，实现非刚性医学图像局部精确配准，本发明能够能很好的模拟非刚性图像的局部形变，从而提高配准的精确度及速度、增强配准的鲁棒性，为医学图像融合奠定一个良好的基础。 

附图说明

图1为本发明非刚性医学图像配准算法的流程图； 

图2为本发明非刚性医学图像配准算法的三角形单元； 

图3为本发明所用的浮动图像； 

图4为本发明所用的参考图像； 

图5为本发明全局配准图像； 

图6为本发明配准算法的配准结果图； 

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。 

如图1所示，本实施例包含以下步骤： 

第一步：对两幅待配准图像进行全局配准 

首先通过一阶矩寻找浮动图像和参考图像这两幅图像的质心，然后通过二阶中心矩寻找浮动图像和参考图像的主轴和坐标系的夹角，根据两幅图像的质心的位置差作为平移分量，然后计算各自的主轴方向并根据其差值将两幅图像旋转对齐，得到旋转分量，实现两幅待配准图像全局配准。具体如下步： 

通过经典力学计算方法计算图像的质心如公式(1)所示： 

$X_{\mathrm{center}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{p}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i},$ $Y_{\mathrm{center}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i{p}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i}...\left(1\right)$

其中：x_i、y_i代表像素点，p_i代表像素的灰度值，n代表图像总共的像素点个数，X_center、Y_center为图像的重心坐标。所述的图像的主轴通过以下方式获得： 

$\tan 2\mathrm{\α}=\frac{2{\mathrm{\μ}}_{11}}{({\mathrm{\μ}}_{20}-{\mathrm{\μ}}_{02})}...\left(2\right)$

其中： 

x，y代表图像中任一点的坐标，p(x，y)代表该点的灰度值，m取值范围为[0，2]，n取值范围为[0，2]，α为图像主轴和参考坐标系的夹角。所述的平移和旋转是指：根据两幅图像的质心的位置差作为平移分量，然后计算各自的主轴方向并根据其差 值将两幅图像旋转对齐，得到旋转分量，实现图像全局配准。实验证明，主轴质心法具有自动、快速、易实现等特点。 

第二步：对全局配准图像和参考图像建立三角形网格 

首先，将步骤1)中的全局配准图像和参考图像分别进行高斯平滑滤波，消除图像成像过程中生成的噪声；然后分别对平滑处理后的两幅待配准图像利用Harris算子提取全局配准图像和参考图像的角点，其次，利用冒泡算法对所有的角点按照横坐标由小到大进行排序，取其横坐标的最小值x_min和最大值x_max，并记录下来。同理对角点按照纵坐标由小到大进行排序，取其纵坐标的最小值y_min和最大值y_max，记录下来，由横纵坐标的最小和最大值，形成四个二维空间的点ABCD，它们坐标分别为A点(x_min，y_min)，B点(x_max，y_min)，C点(x_max，y_max)，D点(x_min，y_max)，连接四点形成四边形ABCD。再次，将四边形ABCD的两个邻边AB(x_max-x_min)、BC(y_max-y_min，)进行等距离的平分，将形成四边形网格。然后，分别正对角线的连接四边形的对角线，将形成三角形网格。最后，将三角形网格进行编号分别为1、2、......、M，M为大于1的整数。 

其角点的提取过程如下： 

2.1)将所处理的矩形区域窗w向任意方向移动位移(x，y)，得到对应灰度值的改变量为： 

$E(u,v){|}_{(x,y)}=\underset{x,y}{\mathrm{\Σ}}w(x,y){[I(x+u,y+v)-I(x,y)]}^2...\left(3\right)$

$E(u,v){|}_{(x,y)}\≅\underset{x,y}{\mathrm{\Σ}}w(x,y){[u^2\frac{{\∂I}^2}{\∂x}+2\mathrm{uv}\frac{\∂I}{\∂x}\frac{\∂I}{\∂y}+v^2\frac{{\∂I}^2}{\∂y}]}^2={\mathrm{Au}}^2+2\mathrm{Cuv}+{\mathrm{Bv}}^2$

$=\left[\begin{array}{cc}u& v\end{array}\right]M\left[\begin{array}{c}u\\ v\end{array}\right]...\left(4\right)$

其中：E(u，v)|_(x，y)表示在点(x，y)处移动一个(u，v)小窗口所发生的灰度值的变化情况，I(u，v)表示图像像素点(x，y)的灰度值，(u，v)表示图像的移动变量，w(x，y)为高斯平滑因子。I(x+u，y+v)-I(x，y)表示图像灰度差值， 

表示图像横坐标的梯度值， 

表示图像纵坐标的梯度值， $A={\left(\frac{\∂I}{\∂x}\right)}^2\⊗w,$ $B={\left(\frac{\∂I}{\∂y}\right)}^2\⊗w,$ $C=\left(\frac{\∂I}{\∂x}\frac{\∂I}{\∂y}\right)\⊗w,$ $M=\left[\begin{array}{cc}A& C\\ C& B\end{array}\right],$ $w_{u,v}=\exp [-\frac{1}{2}(u^2+v^2)/{\mathrm{\δ}}^2]$ 为高斯窗口，对图像窗口进行高斯平滑，目的是提高抗噪能力。 

2.2)利用水平、竖直差分算子对图像每个像素进行卷积后求得I_x和I_y，进而求得矩阵  $M=\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$ 中四个元素的值：，其中：I_x为 I_y为 

I_xI_y为 

所述的卷积是指：采用水平差分算子 $\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]$ 对图像中的每个点可以用水平算子进行卷积，以及采用垂直差分算子 $\left[\begin{array}{ccc}1& 1& 1\\ 0& 0& 0\\ -1& -1& -1\end{array}\right]$ 对图像中的每个点可以用垂直算子进行卷积。 

2.3)对矩阵M的四个元素进行高斯平滑滤波，得到更新后的矩阵M， $M^{\′}=\left[\begin{array}{cc}{A}^{\′}& C^{\′}\\ C^{\′}& B^{\′}\end{array}\right]$ 以计算图像中每个像素的对应角点量R， 

2.4)对图像进行分块得到若干个m*n大小的图像块，对图像块中的角点量R值进行由大到小的排序，经阈值筛选后得到图像的角点。 

所述的分块是指：按照长度和宽度方向正整数方式进行图像划分。 

所述的阈值筛选是指：取图像块中所有角点量的中值或平均值为阈值，当该角点量大于阈值且为某邻域的局部极大值时，该角点为图像的角点。 

第三步：对全局配准图像和参考图像建立数学模型 

本发明将非刚性医学图像的形变过程看作二维弹性薄板的伸缩过程，当作用于弹性体的内力与外力相等时，物体停止形变。因此我们将非刚性医学图像的形变过程看作有一个外力迫使它发生形变的过程，当作用于物体的外力为零时，图像形变结束。具体步骤如下： 

假设I₁(x，y)表示为浮动图像的灰度函数，I₂(x，y)表示参考图像的灰度函数，[x，y]∈R²表示像素空间坐标，h表示从浮动图像到目标图像的空间变换，则任意一点的空间变换可以由如下的映射关系来表示： 

其中u₁(x，y)，u₂(x，y)均为关于坐标(x，y)的函数，记((x，y)-u(x，y))表示形变后的源图像，u(x，y)为形变位移场。配准的目的是期望形变的源图像和目标图像之间的相似度尽可能达到最大。本文中相似度用弹性系统总势能函数来表示，如公式(5) 

其中，σ^t表示弹性体的应力，ε表示弹性体的应变，∫_Ωσ^tεdΩ表示弹性系统的形变势能，Ω表示二维图像区域，∫_ΩF_iudΩ表示弹性系统的外力势能，F_i表示作用于非刚性医学图像的外力，u表示图像的形变位移。 

当两幅待配准图像的对应解剖点在空间位置上达到一致时，则系统的总势能达到最小。 

其次，我们将作用于待配准图像的外力看作是浮动图像和参考图像的灰度差值，如公式6所示表示两幅图像的灰度差值。 

F_i＝I₁((x，y)+u(x，y))-I₂(x，y)..........................................(6) 

则得到：E_总＝∫_Ωσ^tεdΩ+∫_Ω(I₁((x，y)+u(x，y))-I₂(x，y))udΩ......(7) 

当作用于系统外力为零时系统的总势能达最小，图像形变结束，即图像实现局部精确配准。 

第四步：对全局配准图像和参考图像建立有限元方程 

在对步骤3)进行求解的过程中。由于对整幅图像进行建模，使得方程求解过程比较复杂，所以将浮动图像和参考图像进行离散化，将图像感兴趣的区域划分成规则形的三角形网格，将整幅图像的形变看作为三角形单元形变的总和。如步骤2)所述。将步骤2)全局配准图像和参考图像这两幅待配准图像网格单元建立有限元方程，首先，对两幅待配准图像的网格单元进行匹配，将匹配后的网格单元，进行编号1、2、......、N，N为大于1的整数：然后图像的总势能看作为三角形势能的总和，建立有限元方程。 

根据有限元分析中的变分原理和加权余量法将连续的求解图像区域离散为一组有限且按一定方式相互联结在一起单元的组合体。在每个单元内选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，通过插值函数计算出各个单元内场函数的近似解，从而实现两幅图像的配准。 

如步骤3)所述当系统的总势能达到最小时，非刚性医学图像的配准将结束。根据第二步我们可以对非刚性医学图像建立三角形网格，这样系统的总位等于有限三角单元的势能之和，如公式8所示： 

其中，m表示三角形单元的数目，i表示第几个三角形如图2所示。假设非刚性医学图像三角形单元区域的形变u(x，y)量看作较小且光滑，形变后的灰度信息可以表示为如公式9所 示： 

$I_1((x,y)+u(x,y))\&cong;I_1((x,y)+<\&dtri;I_1((x,y),u(x,y)>...<9>$

由此可知非刚性医学图像总势能如公式10所示： 

第五步：求解有限元方程 

如步骤2)所述我们将非刚性医学图像划分为不同的三角形单元，形成三角形网格。 

三角形单元可以由三角形单元线性积分近似，则三角形单元内任意点位移可以由三角形单元的三个节点位移和三角形单元形函数加权和得到如公式11所示： 

$u_i=N_i^{\mathrm{el}}{u}_i^l+N_i^{\mathrm{em}}{u}_i^m+N_i^{\mathrm{en}}{u}_i^n...\left(11\right)$

其中u_i表示三角形单元内的任意点的位移， 

表示三角形单元的形函数，S为三角形单元的面积，l，m，n为三角形单元的三节点其坐标为l(x^l，y^l)，m(x^m，y^m)，n(xⁿ，yⁿ)，则：a^el＝x^myⁿ-xⁿy^m，b^el＝y^m-yⁿ，c^el＝xⁿ-x^m。公式6中的其余系数可由循环置换角标l，m，n的坐标而得到。 

如步骤4)所述当系统的总势能达到最小时，非刚性医学图像的配准将结束。则有：  ${\&Integral;}_{\mathrm{\&Omega;}}\underset{j=1}{\overset{i=3}{\mathrm{\&Sigma;}}}(B_i^{{\mathrm{el}}^t}{\mathrm{DB}}_j^{\mathrm{el}}+N_i^{\mathrm{el}}{\&dtri;I}_1^t{\&dtri;I}_1{N}_j^{\mathrm{el}})u_j^{\mathrm{el}}\mathrm{d\&Omega;}={\&Integral;}_{\mathrm{\&Omega;}}(I_1-I_2){\&dtri;I}_1{N}_i^{\mathrm{el}}\mathrm{\&Omega;},$ 将非刚性医学图像配准问题看作是 

求解在各个单元节点的最佳位移，依次迭代可求的整个非刚性医学图像的形变，实现非刚性医学图像的配准。 

实验结果分析 

表1.1全局配准参数 

质心X轴的坐标	222.485
		质心Y轴的坐标	254.008
水平偏移量	-0.859382
		垂直偏移量	0.861449
角度偏移量	0.0120708

[0075]  由表1.1可知图像全局配准是水平、垂直偏移量、角度偏移量。 

表1.2几种配准算法配准结果对比表 

由表1.2可知本发明基于有限元模型的非刚性医学图像配准算法在配准的精确度有所提高，克服了以往配准算法精确度不高的缺点。 

Claims (1)

1.一种基于有限元模型非刚性医学图像的配准方法，其特征在于：包含如下步骤：

1)全局配准：首先通过图像的一阶矩寻找浮动图像和参考图像的质心；再次，通过浮动图像和参考图像的二阶中心矩寻找所述浮动图像和参考图像的主轴和坐标系的夹角；最后，通过平移和旋转使得浮动图像和参考图像的质心和主轴对齐，得到旋转分量，实现浮动图像和参考图像的全局配准；

2)构建浮动图像和参考图像的网格：首先，将步骤1)中的全局配准图像和参考图像分别进行高斯平滑滤波，消除图像成像过程中生成的噪声；然后分别对平滑处理后的两幅待配准图像利用Harris算子提取图像的角点，对其角点进行排序，取角点的横坐标中最小值、最大值及纵坐标中最小值、最大值并记录下来，这样就形成四个点，连结四个点形成一个四边形，使得四边形包含所有的角点；等分四边形，使得四边形的两组对边平分成等距离，形成四边形网格，将等分后的四边形按照相同的方向连接对角线形成三角形网格，对三角形网格进行编号1、2、……、M，M为大于1的整数；

3)建立数学模型：将全局配准图像和参考图像的形变过程看作弹性的伸缩过程，利用线性弹性模型中的最小势能约束每个匹配后的形变，外力为全局配准图像和参考图像两幅图像的灰度差值，当外力为零时有限元单元的势能到最小，图像完成精确配准；

4)建立有限元方程：在对步骤3)进行求解的过程中，将全局配准图像和参考图像进行离散化，将图像感兴趣的区域划分成规则形的三角形网格，将整幅图像的形变看作为三角形单元形变的总和，对步骤2)中全局配准图像和参考图像这两幅待配准图像网格单元建立有限元方程，首先，对两幅待配准图像的网格单元进行匹配，将匹配后的网格单元，进行编号1、2、……、N，N为大于1的整数，然后图像的总势能看作为三角形势能的总和，建立有限元方程，

二维图像区域，∫_ΩF_iudΩ表示弹性系统的外力势能，F_i表示作用于非刚性医学图像的外力，u表示图像的形变位移，m表示三角形单元的数目，i表示第i个三角形；

5)求解方程：求解步骤4)中的有限元方程，在每个单元内选择至少一个节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，通过插值函数计算出各个单元内场函数的近似解；利用三角形单元的灰度差值的平方和作为外力，当总势能最小时，外力为零，来求解方程，得到全局配准图像的偏移量，实现图像的局部精确配准。

Images