基于自适应三角形网格的非刚性医学图像配准方法

摘要

一种图像处理技术领域的基于自适应三角形网格的非刚性医学图像配准方法，该方法采用分级配准的思想，首先对待配准图像和参考图像进行全局配准；再次利用图像角点来约束全局配准图像和参考图像感兴趣区域生成不规则的三角形网格，在图像感兴趣区域形变较大的地方生成形状较小、数目较多的三角形单元，图像感兴趣区域形变较小的地方生成形状较大、数目较少的三角形单元生成的三角形网格根据图像内容的变化而变化；最后利用近刚性配准的思想实现图像局部精确配准。本发明能够有效的提高非刚性医学图像配准的精确度及增强配准的抗噪能力。

说明书

技术领域

本发明涉及一种图像处理技术领域的方法，具体是一种基于自适应三角形网格的非刚性医学图像配准方法。

背景技术

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

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

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

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

传统的Harris方法在提取图像的角点时需要人为地设置一定的阈值，作为提取角点的一个约束条件，这样会造成在提高图像的阈值时，提取角点的数目少，降低阈值时提取角点数目多，使得提取的角点分布不均匀及易产生聚簇、冗余现象。该方法不能根据图像本质属性自动选取图像的角点。

经过对现有技术的检索发现，王崴、唐一平、时冰川等2008年10月在光学精密工程期刊上发表了文章《一种改进的Harris角点提取方法》，该文章阐述的方法在角点提取时，需要根据实验者以往的实验经验设置一定的阈值，这样既费时又费力且提取角点分布不均匀，容易遗漏角点或产生伪角点。

发明内容

本发明针对现有技术存在的上述不足，提供一种基于自适应三角形网格的非刚性医学图像配准方法，首先对待配准图像和参考图像进行全局配准；再次利用图像角点来约束全局配准图像和参考图像感兴趣区域生成不规则的三角形网格，在图像感兴趣区域形变较大的地方生成形状较小、数目较多的三角形单元，图像感兴趣区域形变较小的地方生成形状较大、数目较少的三角形单元生成的三角形网格根据图像内容的变化而变化；最后利用近刚性配准的思想实现图像局部精确配准。

本发明是通过以下技术方案实现的，本发明包括以下步骤：

第一步：采用主轴质心法对待配准图像和参考图像进行全局刚性配准：首先通过一阶矩寻找待配准图像和参考图像的质心，然后通过二阶中心矩寻找图像的主轴和坐标系的夹角，再通过平移和旋转使得待配准图像和参考图像的质心和主轴对齐，从而实现图像全局刚性配准。

上述的平移和旋转是指：根据待配准图像和参考图像的质心位置差作为平移分量，计算各自的主轴方向并根据其差值将待配准图像和参考图像旋转对齐，得到旋转分量，实现图像配准。

第二步：采用改进的Harris方法提取图像的角点。

第三步aarr：根据Hausdorff角点集匹配方法对全局配准图像和参考图像对应的角点进行匹配。

第四步：根据Delaunay性质来约束角点集，将图像感兴趣区域划分为不规则的三角形网格。

第五步：对待配准图像和参考图像的三角形网格进行匹配排序并标号，利用仿射变换对两个三角形网格集合进行形变，利用相关系数作为相似性度量。

本发明抗噪声能力强，能够很好处理一些随机信号的干扰。

附图说明

图1为本发明流程示意图。

图2为本发明不规则三角形网格单元示意图。

图3表示实施例中形成三角形网格示意图。

图4为实施例效果示意图。

其中：(a)为待配准图像、(b)为参考图像、(c)为基于金字塔配准方法效果图、(d)为实施例效果图、(e)为含噪声的待配准图像、(f)为含噪声的参考图像、(g)为基于金字塔配准方法效果图、(h)为实施例效果图。

具体实施方式

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

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

第一步：采用主轴质心法对待配准图像和参考图像进行全局刚性配准：首先通过一阶矩寻找待配准图像和参考图像的质心，然后通过二阶中心矩寻找图像的主轴和坐标系的夹角，再通过平移和旋转使得待配准图像和参考图像的质心和主轴对齐，从而实现图像全局刚性配准。

上述的平移和旋转是指：根据待配准图像和参考图像的质心位置差作为平移分量，计算各自的主轴方向并根据其差值将待配准图像和参考图像旋转对齐，得到旋转分量，实现图像配准。

第二步：采用改进的Harris方法提取图像的角点，具体步骤如下：

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

$>{E(u,v)|}_{(x,y)}=\underset{x,y}{\Sigma }w(x,y){[I(x+u,y+v)-I(x,y)]}^2---\left(1\right)$>

$>{E(u,v)|}_{(x,y)}\cong \underset{x,y}{\Sigma }w(x,y){[u^2\frac{{\partial I}^2}{\partial x}+2\mathrm{uv}\frac{\partial I}{\partial x}\frac{\partial I}{\partial y}+v^2\frac{{\partial I}^2}{\partial 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(2\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)表示图像灰度差值，表示图像横坐标的梯度值，表示图像纵坐标的梯度值，为高斯窗口，对图像窗口进行高斯平滑，目的是提高抗噪能力。

2.2)利用水平、竖直差分算子对图像每个像素进行卷积后求得I_x和I_y，进而求得矩阵中四个元素的值：，其中：I_x为I_y为I_xI_y为所述的卷积是指：采用水平差分算子对图像中的每个点用水平算子进行卷积，

以及采用垂直差分算子对图像中的每个点用垂直算子进行卷积。

2.3)对矩阵M的四个元素进行高斯平滑滤波，得到更新后的矩阵M，以计算图像中每个像素的对应角点量R，

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

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

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

第三步：根据Hausdorff角点集匹配方法对待配准图像和参考图像对应的角点进行角点匹配，具体步骤包括：

3.1)待配准图像和参考图像对应的角点集S₁和S₂之间的Hausdorff距离为H(S₁，S₂)＝max(h(S₁，S₂)，h(S₂，S₁))，其中：全局配准图像的角点有限集合为S₁＝{a₁ a₂ a₃.....a_p-2 a_p-1 a_p}，参考图像的角点有限集合为S₂＝{b₁ b₂ b₃....b_q-2 b_q-1 b_q}，h(S₁，S₂)为从集合S₁到集合S₂的单向Hausdorff距离，即前向Hausdorff距离，h(S₂，S₁)为从集合S₂到集合S₁的单向Hausdorff距离，即后向Hausdorff距离且||·||是点与点之间的欧式距离，其中h(S₁，S₂)＝d，则表示S₁中所有点到S₂中点的距离不超过d。

3.2)利用Hausdorff距离作为匹配准则，通过穷尽搜索的方法，在图像中的所有位置上移动模板，并求取模板与对应图像中被匹配区域上点集的Hausdorff距离。

3.3)以点集的Hausdorff距离的平均值作为度量t，并删除Hausdorff距离大于度量t的角点，则余下的角点对即为所要搜寻的角点对，具体公式如下：其中：N(S₁)为删除出格点后的累加点数，H(S₁，S₂)＝max(h(S₁，S₂)，h(S₂，S₁))，H(S₁，S₂)中的最小值即为正确的匹配角点对。

第四步：如图3所示，根据Delaunay性质来约束角点集，采用逐点插入法将图像划分为不规则的三角形网格，具体步骤如下：

4.1)根据匹配后的角点坐标值确定一个能包含所有角点的矩形，将矩形分成包含两个三角形的三角形网络，该两个三角形分别作为初始Delaunay三角形；

例如，将参考图像中的匹配角点按照横坐标从小到大进行排序，取其最小值x_min、最大值x_max，同理，寻找纵坐标的最小值y_min，最大值y_max，由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)。Delaunay性质将矩形分为两个Delaunay三角形(如图1左图所示)。

4.2)将角点集中一个未处理点插入三角形网络内，其中(x_i，y_j)为角点p的坐标；

4.3)在三角形网络中找出包含的三角形，把与这个三角形的三个顶点相连，生成三个新的三角形；

4.4)应用Delaunay性质中的最大最小角和每个三角形的外接圆不能包含角点集中的其它点性质，对三角形网络进行更新；

4.5)重复步骤4.2)-步骤4.4)，直到所有的角点都被插入三角形网络中；

4.6)最后删掉包含矩形步骤4.1)中矩形的四个顶点的三角形并得到图像的三角形网格。

第五步：对全局配准图像和参考图像的三角形网格进行匹配排序并标号，利用仿射变化对两个三角形网格集合D1和D2进行形变并利用相关系数测试相似程度，具体步骤包括：

5.1)利用最小二乘法估计全局配准图像和参考图像的三角形网格中的每对三角形单元中具有一对一的匹配对应关系的三对顶点，即三对匹配点的仿射变换参数，即使一个点集中的顶点经过该变换后的坐标与另一个点集中对应点的坐标的欧式距离的平方和S(tx，ty，s，θ)最小变换，其中：令r_i为p_i的坐标与q_i变换后的坐标Z(q_i)的差：

$>r_i=\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)+\left(\begin{array}{cc}s\cos \theta & -s\sin \theta \\ s\sin \theta & s\cos \theta \end{array}\right)\left(\begin{array}{c}{x}_{\mathrm{qi}}\\ y_{\mathrm{qi}}\end{array}\right)-\left(\begin{array}{c}{x}_{\mathrm{pi}}\\ y_{\mathrm{pi}}\end{array}\right)$>

$>=\left(\begin{array}{cccc}1& 0& x_{\mathrm{qi}}& -y_{\mathrm{qi}}\\ 0& 1& y_{\mathrm{qi}}& x_{\mathrm{qi}}\end{array}\right)\left(\begin{array}{c}{t}_x\\ t_y\\ s\cos \theta \\ s\sin \theta \end{array}\right)-\left(\begin{array}{c}{x}_{\mathrm{pi}}\\ y_{\mathrm{pi}}\end{array}\right)$>

其中：w＝(tx ty s cosθ ssinθ)^T，

则上述三对匹配点对的坐标差表示为：

$>R=\left(\begin{array}{c}{r}_1\\ r_2\\ r_3\end{array}\right)=\left(\begin{array}{c}{D}_{q1}w-p_1\\ D_{q2}w-p_2\\ D_{q3}w-p_3\end{array}\right)=\left(\begin{array}{c}{D}_{q1}\\ D_{q2}\\ D_{q3}\end{array}\right)w-\left(\begin{array}{c}{p}_1\\ p_2\\ p_3\end{array}\right)=D·w-P$>

所以将S(tx，ty，s，θ)表达为：

$>S\left(w\right)=\underset{i=1}{\overset{3}{\Sigma }}{r_i}^T{r}_i=\left(\begin{array}{ccc}{r_1}^T& {r_2}^T& {r_3}^T\end{array}\right)\left(\begin{array}{c}{r}_1\\ r_2\\ r_3\end{array}\right)$>

$>=R^{T}R={(\mathrm{Dw}-p)}^T(\mathrm{Dw}-p)$>

$>=w^T{D}^T\mathrm{Dw}-p^T\mathrm{Dw}-w^T{D}^{T}p+p^{T}p$>

$>=w^T{D}^T\mathrm{Dw}-{2p}^T\mathrm{Dw}+p^{T}p$>

S(w)对w求导得：当得到w`：w`＝|D^T d|^-1D^T P；当S(w)对w的二阶导数大于零时，则上式的解S(w)取最小值；求S(w)对w的二阶导数得到：

$>\frac{d^{2}S\left(w\right)}{{\mathrm{dw}}^2}={2D}^{T}D>0;$>

5.2)当全局配准图像和参考图像的相关系数达到最大时，表示全局配准图像和参考图像处于最佳配准位置，利用相关系数作为三角形对的相似性度量：

当全局配准图像的第一个三角形和参考图像的第一个三角形所含信息的重叠程度达到最大时则三角形对相似度量最大，实现图像局部配准；当相似性测度满足一定的误差ε时，则进行下一个匹配三角形对的形变，否则继续进行参数优化，直到相似性测度满足一定的误差ε时，进行下一对三角形形变。以此类推，当第N对匹配三角形进行形变满足一定的误差ε时，且执行到最后一对匹配三角形，则配准结束。若不是执行到最后一对三角形则继续上一次的操作，直到整幅图像的每对匹配三角形执行完为止。

所述的相关系数为：其中：f_i，j、g_i，j分别为全局配准图像和参考图像上三角形的灰度值，分别为两幅图像中三角形的平均灰度值。

实验结果分析：

1)配准的精度、准确性：先将一幅医学图像进行不同角度的翻转、拉伸、平移，然后分别用一种非刚性配准方法(边缘金字塔方法)和本发明的方法进行比较；验证参数包括：灰度差值平方和测度、相关系数、归一化互信息。

如图4(a)、(b)、(c)、(d)及表1所示，本发明配准的效果为比较好、精确度高。

表1实验结果对比表

  配准方法  灰度差值平方和测度  相关系数  归一化互信息  基于边缘金字塔  55.7536  0.894683  1.637733  基于自适应三角网格  32.4758  0.970746  1.874925

2)配准的抗噪性：先将一组待配准的图像加入高斯白噪声，然后用边缘金字塔方法和本发明的方法对加入高斯白噪声的图像进行配准，比较两种方法的配准效果。

如图4中(e)、(f)、(g)、(h)和表2所示，本发明抗噪声能力强，能够很好滤除一些随机干扰信号。

表2加噪声实验结果对比表

  配准方法  灰度差值平方和测度  相关系数  归一化互信息  基于边缘金字塔  65.0864  0.851952  1.533543  基于自适应三角网格  48.6307  0.941825  1.845935