一种基于图像边缘的多分辨非刚性头部医学图像配准方法

摘要

一种基于图像边缘的多分辨非刚性头部医学图像配准方法属于信息融合领域。本发明在配准过程中利用由粗到精的搜索变换参数。首先由小波变换检测图像边缘，在边缘图像基础上构造边缘金字塔图像，然后在边缘金字塔图像的最小尺度层(即最粗糙层)进行较大范围的搜索，寻找配准两幅图像的最佳变换，以后每一层以前一层的搜索结果为中心，缩小搜索范围继续搜索，直到最大尺度层。这种方法与直接用原始图像进行配准相比，可以减少计算量，提高准确度并能够配准不同空间分辨率的图像。

说明书

技术领域

本发明涉及一种基于小波变换的头部医学图像配准方法，是信息融合领域中一项多分辨图像配准方法，在远程医疗系统、医学图像处理等系统中均可有广泛应用。

背景技术

图像配准是图像分析中的一项非常重要的技术，主要通过寻找某种变换，使两幅图像的对应点达到空间位置的一致，在医学诊断过程中，由于存在不同模式图像表现不同性质的物理机制、患者摆位的差异、成像参数的变化、不同成像设备间空间分辨率不同等现实问题，所以单单凭借医生手动将两张或两组不同模式的图像在空间上做对准受到很多局限，且常带有较大的主观性，其可靠性往往不高，不可避免地会产生误差。特别是在定向放射外科和颅脑手术可视化等应用领域，对图像配准的精度要求很高，使得医学图像配准成为一项必要而又相当困难的任务。

图像配准技术的基本过程主要分为三个步骤：

①寻找图像中的对应特征量并提取出来；

②根据特征量寻找最佳匹配变换；

③利用得到的最佳变换对待配准图像进行变换和参考图像进行匹配。其中前两步是配准过程的关键，也是配准算法研究的核心内容。

发明内容

本发明的目的在于提供一种适合于灰度图像特点的配准算法，加强配准算法的针对性，以提高配准后的图像质量，达到理想的实用效果。

在配准过程中利用由粗到精的搜索变换参数。首先由小波变换检测图像边缘，在边缘图像基础上构造边缘金字塔图像，然后在边缘金字塔图像的最小尺度层(即最粗糙层)进行较大范围的搜索，寻找配准两幅图像的最佳变换，以后每一层以前一层的搜索结果为中心，缩小搜索范围继续搜索，直到最大尺度层。这种方法与直接用原始图像进行配准相比，可以减少计算量，提高准确度并能够配准不同空间分辨率的图像。全像素分层匹配时，往往出现对应的同名点不止一个问题，难以实现全局点的匹配。本方法采用基于小波分解的边缘图像金字塔，利用了互信息作为相似性度量准则，分层实现整幅图像的配准。

本发明技术方案的创新点在图像融合之前的预处理和配准规则，为实现这样的目的，本发明的技术方案是：基于小波分解的边缘图配准算法，其特征在于：

(1)使用小波分解进行边缘图像金字塔构建，小波分解是一种无损变换，并且它的直流分量是光滑的，每一层都由一个直流分量和三个特征分量组成，可以利用特征分量进行特征匹配。

(2)预处理后图像每层上的数据与原始图像大小一样，且可以重建，所以在内存中可以释放原始图像，这样，金字塔分层后所需要的内存大小没有增加。与基于互信息的配准方法比较，匹配更加准确，整幅图像的每一层低分辨率的划分，由于像素点仍然相对较多，容易出现匹配结果可能不止一个位置的情况，而边缘图像金字塔可以很好的解决这个问题，减少匹配误差提高配准精度。

本发明提供的一种基于图像边缘的多分辨非刚性头部医学图像配准方法，其特征在于，包括以下步骤：

1)图像边缘提取

设二维图像信号f(x，y)表示图像的像素值，

设η(x，y)为二维平滑函数并且满足：

$\int \int \eta (x,y)\mathrm{dxdy}=1---\left(1\right)$

$\underset{x,y\to \infty }{\lim }\eta (x,y)=0---\left(2\right)$

对平滑函数η(x，y)分别求x方向和y方向的偏导数，则有：

${\mu }^1(x,y)=\frac{\partial \eta (x,y)}{\partial x};$${\mu }^2(x,y)=\frac{\partial \eta (x,y)}{\partial y}$其中μ¹(x，y)和μ²(x，y)看作二维小波函数。

则f(x，y)对应的小波变换为：W^xf(x，y)＝f(x，y)*μ¹

                           W^yf(x，y)＝f(x，y)*μ²

对图像进行平滑，平滑后的图像g(x，y)为：

g(x，y)＝η(x，y)*f(x，y)    (3)

对平滑后的图像求一阶微分，得：

$\frac{\partial g(x,y)}{\partial x}=\frac{\partial }{\partial x}[\eta (x,y)*f(x,y)]$

$=\frac{\partial g(x,y)}{\partial x}*f(x,y)---\left(4\right)$

$\frac{\partial g(x,y)}{\partial y}=\frac{\partial }{\partial y}[\eta (x,y)*f(x,y)]$

$=\frac{\partial g(x,y)}{\partial y}*f(x,y)---\left(5\right)$

式(4)和(5)的右端实际上是函数f(x，y)的两个小波变换，即

$\frac{\partial g(x,y)}{\partial x}=W^{x}f(x,y)=f(x,y)*{\mu }^1(x,y)---\left(6\right)$

$\frac{\partial g(x,y)}{\partial y}=W^{y}f(x,y)=f(x,y)*{\mu }^2(x,y)---\left(7\right)$

一阶导数的模极大值对应图像的边缘点，这两个一阶导数等于f(x，y)的两个小波变换，这两个小波变换的模的极大值就对应了图像的边缘点，从而

由上述方法，得到边缘图像I(x，y)和I′(x，y)，进行图像配准，其中I(x，y)为参考图像对应的边缘图像，I′(x，y)为待配准图像对应的边缘图像；

2).图像配准

2.1)将上述到边缘图像I(x，y)和I′(x，y)变换为各分辨率层的边缘金字塔图像I_k(x，y)和I′_k(x，y)，

过程如下：通过将原始边缘图像I(x，y)和I′(x，y)相邻的2×2个像素算术平均为一个像素构造第一级边缘金字塔图像I₁(x，y)和I′₁(x，y)，然后在第一级边缘金字塔图像I₁(x，y)和I′₁(x，y)的基础上按上述方法构造第二级图像I₂(x，y)和I′₂(x，y)，依此类推构成所需层的边缘金字塔图像I_k(x，y)和I′_k(x，y)；

2.2)在分解成分辨率由低到高的各层边缘金字塔子图像集I_k(x，y)和I′_k(x，y)后，由低一级的分辨率层I_k(x，y)和I′_k(x，y)的搜索结果作为高一级的分辨率层I_k-1(x，y)和I′_k-1(x，y)的初始值，依次类推实现原始图像的配准；

在第k层通过方向加速法进行搜索，设置搜索步长，将参考图像和待配准图像的低频边缘图像进行配准运算，得与I_k(x，y)上点P_K(x，y)在I′_k(x，y)上的最佳对应点P′_k(x，y)，确定图像间的变换参数，并作为下一层搜索的粗略位置；

其中k为边缘金字塔图像分解的第k层，P_K(x，y)为图像I_k(x，y)上的一像素点，x，y为其横纵坐标，P′_k(x，y)为图像I′_k(x，y)上的一像素点；

2.3)在第k-1层配准中以第k层结果P′_k(x，y)作为初始位置，在P′_k(x，y)的邻域内进行配准搜索，在待配准图像中寻找P′_k-1(x，y)使之与参考图像中点P_k-1(x，y)对应，从而得到本层的最佳对应点P′_k-1(x，y)；

2.4)令k＝k-1，重复进行步骤2.3)的配准搜索；

2.5)缩小搜索步长，重复步骤2.3)和步骤2.4，实现逐层搜索直到k＝0，实现k＝0即最顶层上，得到的原始图像上P₀(x，y)和其最佳对应点P′₀(x，y)，即为的最终搜索结果，完成配准过程。

本发明的图像融合方法具有如下有益效果：

(1)本发明算法与一般的基于边缘检测配准方法比较，减少了搜索位置，从而大大提高了配准速度，对于一幅n×n个像素的图像，采用全像素边缘配准算法时，理论上搜索需要的时间为O(n)次，而本发明中算法为O(log n)次(0(x)随x的值单调变化)，扣除预处理构建边缘金字塔图像所花费的时间，仍减少很多时间。

(2)本发明算法与基于全像素金字塔配准方法比较，配准更为准确，全像素金字塔配准构建各层图像时，由于整幅图像中像素点仍然相对较多，往往出现同名像素点不止一个位置的情况，容易造成死循环，而本发明中算法可以很好的解决这个问题，减少了配准误差，提高了配准精度。

(3)本发明使用小波分解进行图像边缘金字塔的构建，小波变换是一种无损变换，并且它的直流分量是光滑的，图像边缘金字塔的每一层都由一个直流分量和三个特征分量组成，可以利用特征分量进行特征匹配。图像金字塔每层上的数据与原始图像大小一样，且可以重建，所以在内存中可以释放原始图像，这样边缘金字塔图像分层后所需要的内存大小没有增加，同时小波变换具有快速算法，配准计算量小速度快。

附图说明

图1本发明实施例子流程图。

图2图像结构金字塔。图中显示的是对原始图像进行的二级金字塔分解示意图。

图3图像配准算法流程图。

图4图像配准效果图。图4.1是CT1109参考图像，图4.2是PD1110待配准图，图4.1和图4.2均是Visual Human图像库中的标准图像。图4.3是将图4.1和图4.2利用本发明配准后的图像。

图5是矢量“共轭程度”示意图。

具体实施方式

图像中各种结构的边缘特征的研究对于模式识别是最重要的。边缘是图像最基本的特征，图像的主要信息大部分都存在于图像的边缘中，这同人们从一幅画的轮廓辨别物体相似。边缘检测在计算机视觉、图像分析等应用中起着重要的作用，是图像分析与识别的重要环节。这是因为图像的边缘包含了用于识别的有用信息，它为人们描述或识别目标和解释图像提供了重要的特征参数。因此，我们边缘检测是图像分析和模式识别的主要特征提取手段。

1.图像边缘提取

小波变换是检测突变信号的强有力工具，能够刻画各种不同频率分量的信号，基于小波变换的边缘检测就是将原始图像进行小波变换，将其分解在不同频段，找出高频部分模的极大值后，进行筛选从而得到图像边缘。

基于小波变换的边缘检测就是用一平滑函数η(x，y)，在不同的尺度下平滑所检测的信号，根据一次或二次微分找出它的突变点(一次微分的极大值点和对应的二次微分的零交叉点和平滑后信号的拐点)，其中所选择的小波函数$\mu \left(t\right)=\frac{\mathrm{d\eta }\left(t\right)}{\mathrm{dt}},$根据小波变换系数极值进行边缘检测。

设二维图像信号f(x，y)，g(x，y)，其中f(x，y)、g(x，y)分别表示图像的像素值，

设η(x，y)为二维平滑函数并且满足：

$\int \int \eta (x,y)\mathrm{dxdy}=1---\left(1\right)$

$\underset{x,y\to \infty }{\lim }\eta (x,y)=0---\left(2\right)$

对平滑函数η(x，y)分别求x方向和y方向的偏导数，则有：

${\mu }^1(x,y)=\frac{\partial \eta (x,y)}{\partial x};$${\mu }^2(x,y)=\frac{\partial \eta (x,y)}{\partial v}$其中μ¹(x，y)和μ²(x，y)可以看作二维小波函数。

对应的小波变换为：W^xf(x，y)＝f(x，y)*μ¹

                  W^yf(x，y)＝f(x，y)*μ²

对图像进行平滑，平滑后的图像为：

g(x，y)＝η(x，y)*f(x，y)    (3)

对平滑后的图像求一阶微分，得：

$\frac{\partial g(x,y)}{\partial x}=\frac{\partial }{\partial x}[\eta (x,y)*f(x,y)]$

$=\frac{\partial g(x,y)}{\partial x}*f(x,y)---\left(4\right)$

$\frac{\partial g(x,y)}{\partial y}=\frac{\partial }{\partial y}[\eta (x,y)*f(x,y)]$

$=\frac{\partial g(x,y)}{\partial y}*f(x,y)---\left(5\right)$

式(4)和(5)的右端实际上是函数f(x，y)的两个小波变换，即

$\frac{\partial g(x,y)}{\partial x}=W^{x}f(x,y)=f(x,y)*{\mu }^1(x,y)---\left(6\right)$

$\frac{\partial g(x,y)}{\partial y}=W^{y}f(x,y)=f(x,y)*{\mu }^2(x,y)---\left(7\right)$

一阶导数的模极大值对应图像的边缘点，这两个一阶导数等于f(x，y)的两个小波变换，因此这两个小波变换的模的极大值就对应了图像的边缘点。

由上述方法，得到边缘图像I(x，y)和I′(x，y)，进行图像配准。(I(x，y)为参考图像对应的边缘图像，I′(x，y)为待配准图像对应的边缘图像)

2.图像配准

基于正交小波变换的多分辨率分解能将金字塔图像中的每一层图像都分解为新的低频平滑图像和高频细节图像，其中低频平滑图像集中了原始图像的大部分能量，反映了图像的绝大部分信息，所以我们可以使用不同层金字塔图像中的低频图像进行图像的分层配准。

在图像配准中算法需要具有一定的精度，所以参考图像和待配准图像特征点要严格对应，即要找到参考图像和待配准图像上对应的特征点(同名点)，现行算法多通过对搜索区域的每一个点进行处理，来寻找同名点，而在整体工作中而除了对配准有作用的同名点外，对其他所有像素点的工作可以被视为“无为”工作，既浪费了时间又影响了精度。所以在发明的算法采用加速算法，利用分层搜索策略解决此问题。

本发明所用算法大致分为以下两步：

(1)把参考图像和待配准图像进行小波变换，得到边缘图像I(x，y)和I′(x，y)，将边缘图像变换为各分辨率层的边缘金字塔图像I_k(x，y)和I′_k(x，y)，过程简述如下：通过将原始边缘图像I(x，y)和I′(x，y)相邻的″2×2″个像素平均为一个像素构造第一级边缘金字塔图像I₁(x，y)和I′₁(x，y)，然后在第一级边缘金字塔图像I₁(x，y)和I′₁(x，y)的基础上按上述方法构造第二级图像I₂(x，y)和I′₂(x，y)，依此类推构成所需层的边缘金字塔图像I_k(x，y)和I′_k(x，y)。(其中像素平均为算术平均)

(2)在分解成分辨率由低到高的各层边缘金字塔子图像集I_k(x，y)和I′_k(x，y)后，利用由粗到精的图像搜索方案实现待配准图像和参考图像的配准，也即由低一级的分辨率层I_k(x，y)和I′_k(x，y)的搜索结果作为高一级的分辨率层I_k-1(x，y)和I′_k-1(x，y)的初始值，依次类推实现原始图像的配准。首先在待配准图像的第k层边缘图像中搜索其与参考图像中某特征点的大致位置P′_k(x，y)，在确定P′_k(x，y)后，再以P′_k(x，y)为中心的邻域内搜索上一层对应点P′_k-1(x，y)，依次类推实现原始图像间特征点的对应。

以Visible Human中的部分图像说明本发明使用方法。

图4.1是CT1109参考图像，图4.2是PD1110待配准图，图4.3是用本发明配准后的图像。

图像配准过程可对照图3，具体步骤是：

(1)首先将参考图像图4.1及待配准图像图4.2进行小波分解，分解层次深度为k层，得到各层上的边缘图像，并建立边缘图像金字塔(见图2)。

图像金字塔的构建可利用下面的步骤：

对一个由n×n个像素组成的图像，我们定义最上面的第零层(k＝0即原始图像)为顶层，最下面一层为底层，金字塔图像下一层像素的值是由上一层像素值计算得到的，计算公式如下：

$I_{k+1}[\left(\frac{x+1}{2}\right),\left(\frac{y+1}{2}\right)]=\left[\frac{I_k(x,y)+I_k(x,y+1)+I_i(x+1,y)+I_k(x+1,y+1)}{4}\right]$

其中x，y＝1，2，3，…，2^i-1，[t]表示对t＝0.5进行取整。I_k(x，y)表示第K层边缘金字塔图像中坐标为(x，y)的像素点的灰度值。

(2)在第k层通过改进的方向加速法(Powell)进行搜索，设置搜索步长，依据互信息最大的相似性准则，将参考图像和待配准图像的低频边缘图像进行配准运算，得与I_k(x，y)上点P_K(x，y)在I′_k(x，y)上的最佳对应点P′_k(x，y)，确定图像间的变换参数，并作为下一层搜索的粗略位置。(k为边缘金字塔图像分解的第K层，P_K(x，y)为图像I_k(x，y)上的一像素点，x，y为其横纵坐标，P′_k(x，y)为图像I′_k(x，y)上的一像素点)。

基本Powell法的搜索过程如下：

①假定初始点P_K(x，y)，初始方向矢量组d₁，d₂，...，d_n，误差ε＞0。(ε为设置的个小正数)

②从P_K(x，y)出发，沿d_n，搜索z_k。

③从z_k出发，依次沿d₁，d₂，...，d_n搜索得X₁^(k)X₂^(k)，...，X_n^(k)，其中

${X_j}^{\left(k\right)}={X_{j-1}}^{\left(k\right)}-{\lambda }_j{d}_j(j=1,...,n,{X_0}^{\left(k\right)}=z_k)$

④判断$\left|\right|X_n^{\left(k\right)}-P_k(x,y)\left|\right|\le ϵ?$如果是，$P_k^{\prime }(x,y)\Leftarrow X_n^{\left(k\right)},$停止搜索；否则产生新方向$d=X_n^{\left(k\right)}-Z_k.$

⑤构造新的方向矢量组$d_1\Leftarrow d_2,d_2\Leftarrow d_3,...,d_{n-1}\Leftarrow d_n,d_n\Leftarrow d,$并设置新的初始点$X_0\Leftarrow X_n^{\left(k\right)},$$k\Leftarrow k+1.$返回②

基本Powell法存在一个缺陷，新的方向组有可能是线性相关的，它们不能张成n维空间Rⁿ，以后的搜索都在一个仿射子集进行，如果极小点不落在这个仿射子集中，将永远达不到极小点X⁺。

也就是说，基本Powell按方向替换规则$d_1\Leftarrow d_2,d_2\Leftarrow d_3,...,d_{n-1}\Leftarrow d_n,d_n\Leftarrow d$

产生的新的方向组，有可能导致新的方向组线性相关。实际上没有理由每次只替换第一个方向矢量。改进的Powell法的要点是：用新的方向矢量d替换原来方向组的某一个方向矢量后，使得新的方向组“共扼程度”尽可能的提高，至少不降低。如果替换原来方向组的任何一个方向都不能提高其“共扼程度”，那就不替换，下一轮循环仍用原来的方向组。

矢量之间要么共扼，要么不共扼，方向组矢量的“共扼程度”或“正交程度”可用图5解释，a1与a2显然比b1与b2更接近正交，共扼也是如此。

定理：设A为n阶正定对称矩阵，方向矢量组D＝(d₁，d₂，...，d_n)满足：

$d_j^{T}A{d}_j=1,(j=\mathrm{1,2},...,n)$

则有：

$\left(1\right),\left|\det \left(D\right)\right|\le \sqrt{\frac{1}{\det \left(A\right)}}$

$\left(2\right),\left|\det \left(D\right)\right|=\sqrt{\frac{1}{\det \left(A\right)}}$的充要条件是d₁，d₂，...，d_n关于A两两共轭。

这一定理说明，如果方向向量组D＝(d₁，d₂，...，d_n)满足$d_j^{T}A{d}_j=1,$当|det(D)|达到最大值时，d₁，d₂，...，d_n关于A两两共扼，因此，在改进的Powell算法中，就以此为标准决定新方向d改替换掉原来方向组d₁，d₂，...，d_n中的哪一个矢量，或者不替换，重新搜索。

(3)在第k-1层配准中以第k层结果P′_k(x，y)作为初始位置，，在P′_k(x，y)的8领域内进行配准搜索，在待配准图像中寻找P′_k-1(x，y)使之与参考图像中点P_k-1(x，y)对应，从而得到本层的最佳对应点P′_k-1(x，y)。

(4)令k＝k-1，重复进行步骤(3)的配准搜索，缩小搜索步长，重复步骤(3)和步骤(4)，实现逐层搜索直到k＝0，实现最顶层I(x，y)和I′(x，y)(即原始参考图像和待配准图像)的特征点配准；

(5)在k＝0层上，得到的原始图像上P₀(x，y)和其最佳对应点P′₀(x，y)，即为的最终搜索结果，将待配准图像进行相应的变换，完成配准过程。

表1图像配准效果比较