CN102722890A - 基于光流场模型的非刚性心脏图像分级配准方法 - Google Patents

Abstract

一种图像处理技术领域的基于光流场模型的非刚性心脏图像分级配准方法，通过两幅图像的尺度不变特征向量得到仿射变换系数，并通过仿射变换得到粗配准图像；再利用光流场方法得到粗配准图像的偏移变换后通过插值得到精配准图像。本发明将SIFT特征方法和光流场方法互为补充，前者能为后者提高收敛速度做准备，后者使得配准结果更为精确；实现更好地保持心脏图像特征细节，具有更强的抗噪能力和鲁棒性，得到较为准确的配准结果。本发明所采用的差值方法是将线性差值和中心差分相结合，同时方法采用多分辨率策略实现最终配准。

Description

基于光流场模型的非刚性心脏图像分级配准方法

技术领域

本发明涉及一种图像处理技术领域的方法，具体是一种基于光流场模型的非刚性心脏图像分级配准方法。

背景技术

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

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

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

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

1980年Moravec等首次提出采用角点检测算子来实现立体视觉匹配，在此基础上Harris等对Moravec算子进行改进，Harris角点检测算子具有旋转不变以及缩放不变等许多优良性能，因此广泛应用在各种图像匹配算法中。但它对尺度、视角、照明变化比较敏感，而且抗噪声能力差。2004年，哥伦比亚大学的Lowe提出了一种新的点特征提取算法——尺度不变特征转换算法(S1FT，Scale-invariant feature transform)，该算法匹配能力较强，能提取稳定的特征，可以对旋转、尺度缩放、亮度变化保持不变性，对视角变化、放射变化、噪声保持一定程度的稳定性，成功应用于图像匹配领域。

Horn等人创造性地将二维速度场与灰度相联系，引入光流约束方程，得到光流计算的基本算法。由于配准所求解的位移场与光流场模型所求解的速度场具有相似性，Palos等人将光流场模型引入到了图像配准中，但这些方法都基于Horn模型。Jean-Philipe Thirion提出的“demons-base”算法是一种简单的基于图像灰度信息的弹性配准方法，和19世纪Maxwell的实验原理很相似，即通过判断出待配准图像上各个像素点的运动方向，然后对每个像素点移动来实现弹性配准。

经过对现有技术的检索发现，白小晶、孙怀江、王平安等2007年12月在江南大学学报上发表的文章《基于光流场模型的医学图像配准》，但是该方法直接利用图像的灰度信息，没有充分利用图像局部特征，在一定程度上限制了配准的精度。

中国专利文献号CN101536919A，公开日2009-9-23，记载了一种“心肌声学造影(MCE)图像定量分析的方法”，该技术步骤为：确定光流计算中的能量方程；为MCE图像序列的各帧构造其金字塔表示；采用由粗到精的搜索策略计算图像序列的光流场；选择某心肌节段的ROI区域；确定该心肌节段所选择ROI区域在MCE图像序列每一帧中的位移；确定该心肌节段所选择ROI区域在MCE图像序列每一帧中的位置，并测量其取样点DB值；确定MCE定量分析数学模型中的未知参数；绘制MCE时间-强度曲线。但该技术主要工作集中在对图像的定量分析上，并没有具体涉及到对非刚性对象心脏的配准。

王安娜、薛嗣麟、俞跃、孙静在“基于改进光流场模型的医学图像配准方法”(第15卷第2期2010年2月《中国图象图形学报》)中提出了对原始光流场模型的正则项进行改进，同时引入运动模糊图像复原方法，改进的方法改善了原始光流场模型造成的图像模糊。实验结果表明，基于改进光流场模型的医学图像配准方法配准结果准确，具有较快的配准速度。但该技术直接对原始光流场算法加入正则项改善图像模糊，没有全面地利用图像局部特征，并且没有预先去除噪声的影响。

林存花、陈海峰在“光流场模型用于非刚性医学图像配准”(《电子科技》2011年04期)中提出了采用一种窗口化的光流配准方法，实现了非刚性医学图像的配准。该方法通过调整窗口函数大小来控制对图像造成的模糊程度，同时采用图像的塔式分解，解决了Horn光流场不能处理大位移的问题。但该技术不能保证图像的拓扑结构不发生变化，没有预先去除噪声的影响。

发明内容

本发明针对现有配准技术存在的上述不足，为了弥补以上技术缺陷提出一种基于光流场模型的非刚性心脏图像分级配准方法，具有良好的尺度、旋转、光照不变特性。光流场方法不需要进行特征提取过程，直接利用图像的灰度信息，具有较快的配准速度。两者互为补充，前者能为后者提高收敛速度做准备，后者使得配准结果更为精确。本发明所采用的差值方法是将线性差值和中心差分相结合。同时方法采用多分辨率策略实现最终配准。

本发明是通过以下技术方案实现的，本发明通过两幅图像的尺度不变特征向量得到仿射变换系数，并通过仿射变换得到粗配准图像；再利用光流场方法得到粗配准图像的偏移变换后通过插值得到精配准图像。

本发明具体包括以下步骤：

第一步：提取参考图像和浮动图像的特征点生成SIFT特征向量，并进行特征匹配，具体步骤包括：

1.1)对待配准图像进行归一化处理，进行预滤波以剔除噪声，得到的图像作为高斯图像金字塔的最底层；

1.2)建立尺度空间，进行尺度空间极值检测得到特征点，并对所得到的特征点进行从大到小的排序，取其中间值作为阈值；然后利用Moravec算子计算特征点的响应函数，当响应函数大于阈值，则保留该点为特征点，否则剔除；

1.3)确定特征点的大小和方向并生成归一化的SIFT特征向量。

所述的大小和方向分别是指：特征点的模值以及邻域像素的梯度方向；

所述的归一化是指：将特征向量的长度进行归一化处理。

1.4)根据两幅图像的SIFT特征向量的欧氏距离，按距离比率(distance-ratio)准则进行特征匹配，即对于一幅图像中某一特征点，当另一幅图像中与其最近的欧式距离为d₁，次欧式距离为d₂，且d₁与d₂的比率大于预设阈值时认为匹配成功；

1.5)根据几何距离条件约束消除特征匹配差错值。

所述的根据几何距离条件约束是指：根据几何定理，到空间中一点的距离为定值的点的轨迹是一个圆，到空间中二点距离为定值的点的轨迹是二个点，到不共线三点的距离分别为定值的点是可以唯一确定的。提取的特征点数目较多而且在图像中比较分散，所以这些点不可能同时共线。因此，可以认为一个点t_a(x，y)是特征点s_a(x，y)的正确匹配点的充要条件是t_a(x，y)点到其他特征点的正确对应点t_i(x，y)的距离和特征点s_a(x，y)到其他特征点s_i(x，y)的距离的个数相等。

第二步：通过最小二乘法得到参考图像和浮动图像的仿射变换系数，并由此对浮动图像进行仿射变换，得到粗配准后的图像。

2.1)参考图像和浮动图像的特征点集P和Q中有n对匹配点对{(p_i，q_i)|i＝1，2，……，n且n≥2}，计算仿射变换参数Z(tx，ty，s，θ)，就是使一个点集中的点经过该变换后的坐标与另一个点集中对应的点的坐标的欧式距离的平方和S(tx，ty，s，θ)最小的变换：

$r_i=\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)+\left(\begin{array}{cc}s\cos \mathrm{\θ}& -s\sin \mathrm{\θ}\\ s\sin \mathrm{\θ}& s\cos \mathrm{\θ}\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 \mathrm{\θ}\\ s\sin \mathrm{\θ}\end{array}\right)-\left(\begin{array}{c}{x}_{\mathrm{pi}}\\ y_{\mathrm{pi}}\end{array}\right)=D_{\mathrm{qi}}w-p_i$

其中：r_i为p_i的坐标与q_i变换后的坐标Z(q_i)的差，w＝(tx ty s cosθ s sinθ)^T， $D_{\mathrm{qi}}=\left(\begin{array}{cccc}1& 0& x_{\mathrm{qi}}& -y_{\mathrm{qi}}\\ 0& 1& y_{\mathrm{qi}}& x_{\mathrm{qi}}\end{array}\right);$

则所有匹配点对的坐标差可以表示为：

$R=\left(\begin{array}{c}{r}_1\\ r_2\\ .\\ .\\ .\\ r_n\end{array}\right)=\left(\begin{array}{c}{D}_{q1}w-p_1\\ D_{q2}w-p_2\\ .\\ .\\ .\\ D_{\mathrm{qn}}w-p_n\end{array}\right)=\left(\begin{array}{c}{D}_{q1}\\ D_{q2}\\ .\\ .\\ .\\ D_{\mathrm{qn}}\end{array}\right)w-\left(\begin{array}{c}{p}_1\\ p_2\\ .\\ .\\ .\\ p_n\end{array}\right)=D\·w-P$

将S(tx，ty，s，θ)表达为 $S\left(w\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{r_i}^T{r}_i=\left(\begin{array}{cccccc}{r_1}^T& {r_2}^T& .& .& .& {r_n}^T\end{array}\right)\left(\begin{array}{c}{r}_1\\ r_2\\ .\\ .\\ .\\ r_n\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$

其中：R表示匹配点对的坐标差，T表示矩阵转置，D表示 $\left(\begin{array}{cccc}1& 0& x_{\mathrm{qi}}& -y_{\mathrm{qi}}\\ 0& 1& y_{\mathrm{qi}}& x_{\mathrm{qi}}\end{array}\right),$ w表示(tx ty s cosθ s sinθ)^T，p表示p_i的坐标值。

2.2)求解S(w)最小时的w，记为w`，则有

即w`＝|D^T d|^-1D^TP；进一步使得S(w)对w的二阶导数必须大于零

求解得到w`是使S(w)最小的解：w`＝|D^T D|^-1D^TP，即仿射变换系数；进一步根据得到的变换参数对待配准图像进行仿射变换，得到粗配准图像。

第三步：针对第二步得到粗配准图像和参考图像进行基于光流场方法的迭代；

3.1)将粗配准图像每一位置向量为x_I的像素点I的初始偏移V_I设为0，分别对参考图像和粗配准图像构造级数为K的高斯多分辨率金字塔；

3.2)计算当前金字塔分辨率情况下每一个粗配准图像像素点的位置向量

其中：K为当前分辨率级别，

为通过偏移公式计算得到的偏移变换，N为迭代次数；

3.3)计算两幅图像间的互信息对当前偏移变换

的梯度，当两幅图像间的互信息的增量小于预设阈值时则认为当前分辨率级数下的迭代已经收敛，将分辨率级数K递增1并进入步骤3.4)，否则返回步骤3.2)且迭代次数为N+1，进行下一次迭代，直至两幅图像间的互信息的增量不小于互信息阈值；

3.4)对前一级分辨率级数下得到的偏移变换

进行重采样，得到当前分辨率级数下的初始偏移变换，返回步骤3.2)，并进行当前分辨率级数下的迭代，直至该分辨率下两幅图像间的互信息的增量不小于互信息阈值，并且各个分辨率级数均执行迭代结束。

所述的重采样是指图像变换常运用重采样技术(即插值运算)来近似计算坐标变换后相应像素的灰度值。

第四步：各分辨率级数所对应的迭代结束后，将所得到的偏移变换用于粗配准图像，经插值得到精配准图像。

所述的插值是指：采用线性插值和中心差分方法对粗配准图像进行插值；

本发明中，精配准图像可进一步通过相关性函数进行校核：对插值后的图像以相关系数为相似性测度进行校核，设插值后图像A(x，y)和参考图像B(x，y)，相关性函数为： $C_{A,B}\left(T\right)=\frac{{\mathrm{\Σ}}_{x_1,y_1}^{x_n,y_n}A\left(T(x,y)\right)\·B(x,y)}{\sqrt{{\mathrm{\Σ}}_{x_1,y_1}^{x_n,y_n}{A}^2\left(T(x,y)\right)}\·\sqrt{{\mathrm{\Σ}}_{x_1,y_1}^{x_n,y_n}{B}^2(x,y)}},$ 其中：A(x，y)和B(x，y)代表图像在(x，y)坐标处的灰度值，T代表图像A和图像B坐标变换，当相关系数C_A，B(T)大于等于相关系数阈值则校核认为达到最佳配准，相关系数阈值一般为0.985。若低于该阈值则需要将精配准图像作为浮动图像与参考图像重新运行本方法。

本发明抗噪声能力强，提高了方法的鲁棒性。

附图说明

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

图2为本发明尺度空间极值点的确定示意图。

图3为实施例一效果示意图；

其中：(a)为参考图像、(b)为浮动图像、(c)为参考图像生成特征向量效果图、(d)为浮动图像生成特征向量效果图、(e)为特征匹配效果图、(f)为粗配准后图像、(g)为实施例效果图、(h)为基于特征的配准效果图、(i)为基于B样条的配准效果图。

图4为实施例二效果示意图；

其中：(a)为参考图像、(b)为浮动图像、(c)为参考图像生成特征向量效果图、(d)为浮动图像生成特征向量效果图、(e)为特征匹配效果图、(f)为粗配准后图像、(g)为实施例效果图、(h)为基于特征的配准效果图、(i)为基于B样条的配准效果图。

具体实施方式

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

实施例1

本实施例采用图像大小为256×256的时间序列的心脏CT图像，相隔时间较短，为相邻帧；

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

第一步：首先对待配准图像进行归一化处理，然后将图像扩大至原来的两倍，进行预滤波以剔除噪声，在这里采用中值滤波滤除脉冲干扰和图像扫面噪声。得到的图像作为高斯图像金字塔的最底层，即第1阶的第1层，并生成SIFT特征向量，具体步骤如下：

1.1)建立尺度空间，进行尺度空间极值检测。

利用高斯核函数对待配准图像进行尺度变换，获取图像多尺度下的尺度空间表示序列，然后对这些表示序列进行尺度空间特征提取。

公式(1)所示为二维高斯核，其中σ表示高斯正态分布的方差：

$G(x,y,\mathrm{\σ})=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{e}^{-(x^2+y^2)/2{\mathrm{\σ}}^2}---\left(1\right)$

二维图像I(x，y)在不同尺度空间表示L(x，y，σ)可由图像I(x，y)与高斯核G(x，y，σ)的卷积得到，如公式(2)所示：

L(x，y，σ)＝G(x，y，σ)×I(x，y)                    (2)

其中，L表示尺度空间，(x，y)代表图像I上的点，σ是尺度因子，其值越大表示该图像被平滑得越大；其值越小表示该图像被平滑得越小。选择合适的尺度因子平滑是建立尺度空间的关键。然后将相邻两个尺度的图像相减建立图像的DOG(Difference-of-Gaussian)金字塔。

在这里，高斯金字塔一般选择4阶，每一阶选择5层。高斯金字塔和DOG金字塔的构成如下表1、表2所示：

在上面建立的DOG尺度空间金字塔中，当一个点在DOG金字塔本层以及上下两层的26个近邻像素的值中是最大或最小值，就认为该点是图像在该尺度下的一个特征点，如图2所示。

1.2)定位特征点位置，去除不稳定的边缘相应点和低对比度的特征点。

由于DOG值对噪声和边缘较敏感，因此，在上面DOG尺度空间中检测到局部极值点还要经过进一步的检验才能精确定位为特征点。图像中某些区域可能有较多的像素点变化大于经验值而被作为特征点提出，造成特征点冗余或者特征点簇聚。所以对所得到的特征点进行从大到小的排序，取其中间值赋给阈值。利用Moravec计算特征点的响应函数，当响应函数大于设定的阈值，则该点为特征点。否则不是特征点。

1.3)确定特征点的大小和方向。

利用特征点邻域像素的梯度方向分布特性为每个特征点指定一个或多个方向，使算子具有旋转不变性。特征点的模值和方向公式分别为公式(3)、(4)所示：

$m(x,y)=\sqrt{{(L(x+1,y)-L(x-1,y))}^2+{(L(x,y+1)-L(x,y-1))}^2}---\left(3\right)$

θ(x，y)＝tan^-1((L(x，y+1)-L(x，y-1))/(L(x+1，y)-L(x-1，y))    (4)

其中：(x，y)要确定是哪一阶的哪一层，L为每个特征点各自所在的尺度。在以特征点为中心的邻域窗口内采样，并用梯度方向直方图统计邻域像素的梯度方向。梯度方向直方图的峰值代表了该特征点所处邻域梯度的主方向，即为该特征点的方向。

当存在一个相当于主峰值80％能量的峰值，就把这个方向作为该特征点的辅方向。由此一个特征点可能会被指定具有多个方向，这样可以增强特征匹配的鲁棒性。

所述的邻域窗口是以该特征点为中心，取3×3的窗口。

1.4)生成SIFT特征向量。

将坐标轴旋转为特征点的方向，以特征点为中心取16×16的窗口，分成16个4×4的子块，在每个4×4的子块上计算8个方向的梯度直方图，可以得到一个具有8个方向的种子点。这样一个特征点会产生4×4×8共128维的特征向量。这时候SIFT特征向量已经去除了旋转、尺度变化等几何变形因素的影响，然后对特征向量的长度进行归一化处理，去除光照变化的影响。

第二步：对浮动图像和参考图像对应的特征向量进行特征匹配，具体步骤包括：

2.1)设待匹配的两个图像为M和N，两幅图片经SIFT方法提取特征点后，特征点集合分别为：图像M的特征点集合为

图像N的特征点集合为

其中L_m和L_n分别为图像M和N的特征点数目。使用欧几里得公式计算两个特征点间的距离，当SIFT特征向量的维度为k时，距离为下式：

$d(F_m,F_n)=\sqrt{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{(f_i^{\left(m\right)}-f_i^{\left(n\right)})}^2}---\left(5\right)$

并按距离比率(distance-ratio)准则进行特征匹配，即对于某一特征点，设另一图像中与其最近的欧几里得距离为d₁，次欧几里得距离为d₂，且d₁与d₂的比率为ratio＝d₁/d₂。评判是否与特征点匹配的定义为：

$\left\{\begin{array}{c}\mathrm{ratio}>\mathrm{\ϵ},\mathrm{success};\\ \mathrm{Else},\mathrm{failure}.\end{array}\right.---\left(6\right)$

其中，ε为预先设定的阈值，也就是说当距离比率大于某一阈值时，认为特征点成功匹配；相反，则认为与特征点不匹配。

2.2)根据几何定理消除错配，具体步骤为：通过相似性度量得到潜在匹配对，其中不可避免会产生一些错误匹配，因此需要根据几何限制和其他附加约束消除错误匹配，提高鲁棒性。依据几何方面的知识，到空间中一点的距离为定值的点的轨迹是一个圆，到空间中二点距离为定值的点的轨迹是二个点，到不共线三点的距离分别为定值的点是可以唯一确定的。提取的特征点数目较多而且在图像中比较分散，所以这些点不可能同时共线。因此，可以认为一个点t_a(x，y)是特征点s_a(x，y)的正确匹配点的充要条件是t_a(x，y)点到其他特征点的正确对应点t_i(x，y)的距离和特征点s_a(x，y)到其他特征点s_i(x，y)的距离的个数相等，称该条件为距离条件约束。满足该条件的匹配对被接受，否则被消除。

第三步：粗配准的实现：通过最小二乘法求取两幅待配准图像的相似变换参数，找到浮动图像在参考图像中的位置，对图像进行粗配准，得到粗配准后的图像。

2.1)两幅待配准图像的特征点集P和Q中有n对匹配点对{(p_i，q_i)|i＝1，2，…，n且n≥2}，计算仿射变换参数Z(tx，ty，s，θ)，就是使一个点集中的点经过该变换后的坐标与另一个点集中对应的点的坐标的欧式距离的平方和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 \mathrm{\θ}& -s\sin \mathrm{\θ}\\ s\sin \mathrm{\θ}& s\cos \mathrm{\θ}\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 \mathrm{\θ}\\ s\sin \mathrm{\θ}\end{array}\right)-\left(\begin{array}{c}{x}_{\mathrm{pi}}\\ y_{\mathrm{pi}}\end{array}\right)$

$=D_{\mathrm{qi}}w-p_i$

其中w＝(tx ty s cosθ s sinθ)^T， $D_{\mathrm{qi}}=\left(\begin{array}{cccc}1& 0& x_{\mathrm{qi}}& -y_{\mathrm{qi}}\\ 0& 1& y_{\mathrm{qi}}& x_{\mathrm{qi}}\end{array}\right)$

2.2)n对匹配点对的坐标差可以表示为：

$R=\left(\begin{array}{c}{r}_1\\ r_2\\ .\\ .\\ .\\ r_n\end{array}\right)=\left(\begin{array}{c}{D}_{q1}w-p_1\\ D_{q2}w-p_2\\ .\\ .\\ .\\ D_{\mathrm{qn}}w-p_n\end{array}\right)=\left(\begin{array}{c}{D}_{q1}\\ D_{q2}\\ .\\ .\\ .\\ D_{\mathrm{qn}}\end{array}\right)w-\left(\begin{array}{c}{p}_1\\ p_2\\ .\\ .\\ .\\ p_n\end{array}\right)$

$=D\·w-P$

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

$S\left(w\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{r_i}^T{r}_i=\left(\begin{array}{cccccc}{r_1}^T& {r_2}^T& .& .& .& {r_n}^T\end{array}\right)\left[\begin{array}{c}{r}_1\\ r_2\\ .\\ .\\ .\\ r_n\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$

2.3)求解S(w)最小时的w，即为w`，因此S(w)对w求导可得：

$\frac{\mathrm{dS}\left(w\right)}{\mathrm{dw}}=2D^T\mathrm{Dw}-2D^{T}P$ (式3.8)

要使

可得到w`：

w`＝|D^T d|^-1D^TP                  (式3.9)

当S(w)对w的二阶导数必须大于零，那么上式的解是使得S(w)最小的值。对S(w)对w求二阶导数可以得到：

(式3.10)

所以，w`是使S(w)最小的解：

w`＝|D^T D|^-1D^TP                  (式3.11)

通过对上式进行最小二乘法的求解，得到所要求解的仿射变换系数。通过得到的变换参数对待配准图像进行仿射变换，得到粗配准图像。

第四步：本发明采用的基于光流场模型的Demons配准方法是一种简单的基于图像灰度信息的弹性配准方法，通过判断出待配准图像上各个像素点的运动方向，然后对每个像素点移动来实现弹性配准，具体步骤包括：

假设F是浮动图像，R是参考图像，得到：

V(x)▽(R(x))＝R(x)-R(x)          (7)

F(x)和R(x)分别是图像F和R在坐标x处的灰度值。▽(R(x))是图像R在坐标x处的灰度梯度，V(x)是坐标从图像F到R的偏移。配准的目的就是找到从图像F到R的偏移。进一步可得到：

$V\left(x\right)=\frac{(F\left(x\right)-R\left(x\right))\▿\left(R\left(x\right)\right)}{\left|\right|\▿\left(R\left(x\right)\right)\left|\right|}---\left(8\right)$

但是，当▽(R(x))→0时，上述公式变得很不稳定，导致较大的偏移V，为解决此问题，在分母增加一个分量，得到：

$V\left(x\right)=\frac{(F\left(x\right)-R\left(x\right))\▿\left(R\left(x\right)\right)}{\left|\right|\▿\left(R\left(x\right)\right)\left|\right|+{(F\left(x\right)-R\left(x\right))}^2}---\left(9\right)$

在每一次迭代后，使用高斯滤波来平滑所得到的偏移，使变换正则化，如下式所示：

$V_{n+1}\left(x\right)=G_{\mathrm{\σ}}\⊗(V_n\left(x\right)+\frac{(F\left(x\right)-R\left(x\right))\▿\left(R\left(x\right)\right)}{\left|\right|\▿\left(R\left(x\right)\right)\left|\right|+{(F\left(x\right)-R\left(x\right))}^2})---\left(10\right)$

其中，高斯滤波为高斯滤波中的标准差σ被称为弹性系数，它决定了整个非刚性配准过程。

两幅图像间的互信息对当前偏移变换的梯度，从而使浮动图像朝着使两图像间互信息增大的方向变换，当两图像间的互信息不再增大时，认为它们已经配准。如下式所示：

$V_{n+1}\left(x\right)=G_{\mathrm{\σ}}\⊗(V_n\left(x\right)+\frac{(F\left(x\right)-R\left(x\right))\▿\left(R\left(x\right)\right)}{\left|\right|\▿\left(R\left(x\right)\right)\left|\right|+{(F\left(x\right)-R\left(x\right))}^2}+\mathrm{\α}\▿\mathrm{MI}\left(V_n\left(x\right)\right))---\left(11\right)$

在上式中，▽MI(V_n(x))是两幅图像间互信息在像素点x处对于当前变换的梯度值。α为正常数，表示迭代的步长。

其中，互信息梯度如下：

在Hermosillo的互信息梯度理论中，设浮动图像F以参考图像R为目标进行配准，其局部变换U＝x+V(x)将R中的点x映射成F中的对应点x，将两幅图像间的互信息对当前空间位移向量V的导数定义为互信息的梯度，并通过Parzen窗将图像灰度的联合分布

改写成连续函数。其中I₁为参考图像在x点处的灰度，I₂为浮动图像在x点处的灰度。

变换后的图像F与参考图像R间的互信息定义为

$I\left(V\right)={\∫\∫P}_V^{R,F}(I_1,I_2)\log \frac{P_V^{R,F}(I_1,I_2)}{P^R\left(I_1\right)P_V^F\left(I_2\right)}d{I}_{1}d{I}_2---\left(12\right)$

当变形域V被扰动成V+εH，得到变分形式：

$\frac{\∂I(V+\mathrm{\ϵH})}{\∂\mathrm{\ϵ}}{|}_{\mathrm{\ϵ}=0}=\∫\∫\frac{\∂}{\∂\mathrm{\ϵ}}{\left[P_V^{R,F}(I_1,I_2)\log \frac{P_V^{R,F}(I_1,I_2)}{P^R\left(I_1\right)P_{V+\mathrm{\ϵH}}^F\left(I_2\right)}\right]}_{\mathrm{\ϵ}=0}d{I}_{1}d{I}_2---\left(13\right)$

$=(1+\log \frac{P_V^{R,F}(I_1,I_2)}{P^R\left(I_1\right)P_{V+\mathrm{\ϵH}}^F\left(I_2\right)})\frac{\∂P_{V+\mathrm{\ϵH}}^{R,F}(I_1,I_2)}{\∂\mathrm{\ϵ}}{|}_{\mathrm{\ϵ}=0}d{I}_{1}d{I}_2$

$P=\∫\frac{\∂P_{U+\mathrm{\ϵH}}^F\left(I_2\right)}{\∂\mathrm{\ϵ}}\frac{1}{P_{U+\mathrm{\ϵH}}^F\left(I_2\right)}\left({\∫P}_{U+\mathrm{\ϵH}}^{R,F}(I_1,I_2)d{I}_1\right)d{I}_2---\left(14\right)$

$=\frac{\∂}{\∂\mathrm{\ϵ}}\left[{\∫\∫P}_{U+\mathrm{\ϵH}}^F\left(I_2\right)d{I}_2\right]=0$

图像灰度的联合熵分布由两幅图像重叠的部分来估算，重叠部分的像素数为N。估算过程中使用宽度为δ的2维Parzen窗函数，Parzen窗法是从测量样本直接估计随机变量概率密度的一种精确的非参数估计方法：

$P_{V-\mathrm{\ϵH}}^{R,F}(I_1,I_2)=\frac{1}{N}\underset{\mathrm{\ω}}{\∫}{\mathrm{\Ψ}}_{\mathrm{\δ}}(R\left(x\right)-I_1,F(x+V\left(x\right)+\mathrm{\ϵH}\left(X\right))-I_2)\mathrm{dx}---\left(15\right)$

的导数可以很容易通过计算得到：

$\frac{\∂P_V^{R,F}(I_1,I_2)}{\∂\mathrm{\ϵ}}=\frac{1}{N}\underset{\mathrm{\ω}}{\∫}{\mathrm{\Ψ}}_{\mathrm{\δ}}(R\left(x\right)-I_1,F(x+V\left(x\right)+\mathrm{\ϵH}\left(X\right))-I_2)\▿F(x+V\left(x\right)+\mathrm{\ϵH}\left(x\right))\·H\left(x\right)\mathrm{dx}---\left(16\right)$

让ε＝0，得到：

$\frac{\∂I(V+\mathrm{\ϵH})}{\∂\mathrm{\ϵ}}{|}_{\mathrm{\ϵ}=0}=\∫\∫\∫L_V(I_1,I_2){\mathrm{\Ψ}}_{\mathrm{\δ}}(R\left(x\right)-I_1,F(x+V\left(x\right)+\mathrm{\ϵH}\left(X\right))-I_2)\▿F(x+V\left(x\right))\·H\left(x\right)\mathrm{dxd}{I}_{1}d{I}_2---\left(17\right)$

在公示(17)中，

▽F(x+V(x)表示浮动图像F的灰度梯度。

使用卷积符号

表示，则有：

$\frac{\∂I(V+\mathrm{\ϵH})}{\∂\mathrm{\ϵ}}{|}_{\mathrm{\ϵ}=0}=\frac{1}{N}{\∫}_V[\mathrm{\Ψ}\⊗\frac{\∂L_V}{\∂I_2}](R\left(x\right),F(x+V\left(x\right)))\▿F(x+V\left(x\right))H\left(x\right)\mathrm{dx}---\left(18\right)$

其中，Ψ_δ(I₁，I₂)＝K_δ(I₁)K_δ(I₂)，

$\frac{\∂L_V}{\∂I_2}=\frac{1}{P_U^{R,F}}\frac{\∂P_U^{R,F}(I_1,I_2)}{\∂I_2}-\frac{1}{P_U^F\left(I_2\right)}\frac{\∂P_U^F\left(I_2\right)}{\∂I_2}---\left(19\right)$

由此可以得到互信息梯度，如公示(20)所示。

$\▿I_V=\frac{1}{N}[{\mathrm{\Ψ}}_{\mathrm{\δ}}\⊗\frac{d{L}_V}{d_{I_2}}](R\left(x\right),F(x+V\left(x\right)))\▿F(x+V\left(x\right))---\left(20\right)$

可以认为浮动图像上的像素点x_I只有沿[▽_xI_V(x_I)，▽_yI_V(y_I)]^T方向发生移动时，才能保证两幅图像间的互信息增加。

第五步：精配准的实现，具体步骤如下：

5.1)将浮动图像每一位置向量为x_I的像素点I的初始偏移V_I设为0，或者对图像进行刚性粗配准得到初始偏移。对参考图像和浮动图像构造高斯多分辨率金字塔，级数为K；

5.2)在小于K的当前分辨率级别下，在第N次迭代情况下，对于每一个浮动图像像素点，其位置向量为

偏移

通过公式(11)计算；

5.3)评价

若两幅图像间的互信息的增量小于预设阈值，则认为当前分辨率级下的迭代已经收敛，进入下一分辨率，转入步骤5.4)；否则，进入下一次迭代N+1；

5.4)对前一级分辨率得到的偏移变换，进行重采样，得到当前分辨率级下的初始偏移变换，返回步骤5.2)，并进行当前分辨率级数下的迭代，直至该分辨率下两幅图像间的互信息的增量不小于互信息阈值，并且各个分辨率级数均执行迭代结束。

5.5)各分辨率级数所对应的迭代结束后，将所得到的偏移区域V作用于浮动图像，插值得到配准后的图像，即精配准图像。

上述步骤得到的精配准图像可进一步通过相关性函数进行校核：对插值后的图像以相关系数为相似性测度进行校核，设插值后图像A(x，y)和参考图像B(x，y)，相关性函数为： $C_{A,B}\left(T\right)=\frac{{\mathrm{\Σ}}_{x_1,y_1}^{x_n,y_n}A\left(T(x,y)\right)\·B(x,y)}{\sqrt{{\mathrm{\Σ}}_{x_1,y_1}^{x_n,y_n}{A}^2\left(T(x,y)\right)}\·\sqrt{{\mathrm{\Σ}}_{x_1,y_1}^{x_n,y_n}{B}^2(x,y)}},$ 其中：A(x，y)和B(x，y)代表图像在(x，y)坐标处的灰度值，T代表图像A和图像B坐标变换，当相关系数C_A，B(T)大于等于相关系数阈值则校核认为达到最佳配准，相关系数阈值一般为0.985。若低于该阈值则需要将精配准图像作为浮动图像与参考图像重新运行本方法。

实施例2

本实施例采用图像大小为512×512的时间序列心脏图像，相隔时间相对前组长些，为相隔帧。实验所用的计算机配置是Intel酷睿2双核CPU，内存2G，主频2.4GHz。

实施例2步骤与实施例1相同。

实验结果分析：

1)配准的精度、准确性：先将一幅256*256大小的医学图像进行不同角度的翻转、拉伸、平移，然后分别用两种非刚性配准方法(基于特征方法、基于B样条)和本发明的方法进行比较；验证参数包括：均方根、相关系数、归一化互信息。

如图4(a)、(b)、(c)、(d)及表3所示，本发明配准的效果为比较好、鲁棒性强，现有技术不能有效的模拟非刚性心脏图像的形变过程，配准的误差比较大。

表3实验结果对比表

配准方法	均方根	相关系数	归一化互信息
				基于特征	19.0304	0.98496	1.63253
基于B样条	18.1256	0.98723	1.70825
				本方法	17.033	0.99237	1.81743

2)配准的鲁棒性：将一幅512*512大小的医学图像分别用两种非刚性配准方法(基于特征方法、基于B样条)和本发明的方法进行比较；验证参数包括：均方根、相关系数、归一化互信息。

如图4中(e)、(f)、(g)、(h)和表2所示，得知上述方法可以应用于不同大小的图像，配准效果比现有技术精确，本发明鲁棒性强。

表4实验结果对比表

配准方法	均方根	相关系数	归一化互信息
				基于特征	16.1413	0.97962	1.53356
基于B样条	16.3575	0.98334	1.66754
				本方法	14.7855	0.98903	1.80545

Claims (10)

1.一种基于光流场模型的非刚性心脏图像分级配准方法，其特征在于，通过两幅图像的尺度不变特征向量得到仿射变换系数，并通过仿射变换得到粗配准图像；再利用光流场方法得到粗配准图像的偏移变换后通过插值得到精配准图像。

2.根据权利要求1所述的方法，其特征是，该方法具体包括以下步骤：

第一步：提取参考图像和浮动图像的特征点生成SIFT特征向量，并进行特征匹配；

第二步：通过最小二乘法得到参考图像和浮动图像的仿射变换系数，并由此对浮动图像进行仿射变换，得到粗配准后的图像；

第三步：针对第二步得到粗配准图像和参考图像进行基于光流场方法的迭代；

第四步：各分辨率级数所对应的迭代结束后，将所得到的偏移变换用于粗配准图像，经插值得到精配准图像。

3.根据权利要求1或2所述的方法，其特征是，所述的特征向量通过以下步骤得到：

1.1)对待配准图像进行归一化处理，进行预滤波以剔除噪声，得到的图像作为高斯图像金字塔的最底层；

1.2)建立尺度空间，进行尺度空间极值检测得到特征点，并对所得到的特征点进行从大到小的排序，取其中间值作为阈值；然后利用Moravec算子计算特征点的响应函数，当响应函数大于阈值，则保留该点为特征点，否则剔除；

1.3)确定特征点的大小和方向并生成归一化的SIFT特征向量。

4.根据权利要求3所述的方法，其特征是，所述的大小和方向分别是指：特征点的模值以及邻域像素的梯度方向；所述的归一化是指：将特征向量的长度进行归一化处理。

5.根据权利要求2所述的方法，其特征是，所述的特征匹配是指：根据两幅图像的SIFT特征向量的欧氏距离，按距离比率准则进行特征匹配，即对于一幅图像中某一特征点，当另一幅图像中与其最近的欧式距离为d₁，次欧式距离为d₂，且d₁与d₂的比率大于预设阈值时认为匹配成功；根据几何距离条件约束消除特征匹配差错值。

6.根据权利要求1或2所述的方法，其特征是，所述的仿射变换系数通过以下步骤得到：

2.1)参考图像和浮动图像的特征点集P和Q中有n对匹配点对{(p_i，q_i)|i＝1，2，……，n且n≥2}，计算仿射变换参数Z(tx，ty，s，θ)，就是使一个点集中的点经过该变换后的坐标与另一个点集中对应的点的坐标的欧式距离的平方和S(tx，ty，s，θ)最小的变换：

$r_i=\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)+\left(\begin{array}{cc}s\cos \mathrm{\θ}& -s\sin \mathrm{\θ}\\ s\sin \mathrm{\θ}& s\cos \mathrm{\θ}\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 \mathrm{\θ}\\ s\sin \mathrm{\θ}\end{array}\right)-\left(\begin{array}{c}{x}_{\mathrm{pi}}\\ y_{\mathrm{pi}}\end{array}\right)=D_{\mathrm{qi}}w-p_i$

其中：r_i为p_i的坐标与q_i变换后的坐标Z(q_i)的差，w＝(tx ty s cosθ s sinθ)^T， $D_{\mathrm{qi}}=\left(\begin{array}{cccc}1& 0& x_{\mathrm{qi}}& -y_{\mathrm{qi}}\\ 0& 1& y_{\mathrm{qi}}& x_{\mathrm{qi}}\end{array}\right);$

则所有匹配点对的坐标差可以表示为：

$R=\left(\begin{array}{c}{r}_1\\ r_2\\ .\\ .\\ .\\ r_n\end{array}\right)=\left(\begin{array}{c}{D}_{q1}w-p_1\\ D_{q2}w-p_2\\ .\\ .\\ .\\ D_{\mathrm{qn}}w-p_n\end{array}\right)=\left(\begin{array}{c}{D}_{q1}\\ D_{q2}\\ .\\ .\\ .\\ D_{\mathrm{qn}}\end{array}\right)w-\left(\begin{array}{c}{p}_1\\ p_2\\ .\\ .\\ .\\ p_n\end{array}\right)=D\·w-P$

将S(tx，ty，s，θ)表达为 $\left(w\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{r_i}^T{r}_i=\left(\begin{array}{cccc}{r_1}^T& {r_2}^T& ...& {r_n}^T\end{array}\right)\left(\begin{array}{c}{r}_1\\ r_2\\ .\\ .\\ .\\ r_n\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$

其中：R表示匹配点对的坐标差，T表示矩阵转置，D表示 $\left(\begin{array}{cccc}1& 0& x_{\mathrm{qi}}& -y_{\mathrm{qi}}\\ 0& 1& y_{\mathrm{qi}}& x_{\mathrm{qi}}\end{array}\right),$ w表示(tx ty s cosθ s sinθ)^T，p表示q_i的坐标值；

2.2)求解S(w)最小时的w，记为w`，则有

即w`＝|D<sup>T</sup> d|<sup>-1</sup>D<sup>T</sup>P；进一步使得S(w)对w的二阶导数必须大于零求解得到w`是使S(w)最小的解：w`＝|D^T D|^-1D^TP，即仿射变换系数。

7.根据权利要求1或2所述的方法，其特征是，所述的光流场方法是指：

3.1)将粗配准图像每一位置向量为x_I的像素点I的初始偏移V_I设为0，分别对参考图像和粗配准图像构造级数为K的高斯多分辨率金字塔；

3.2)计算当前金字塔分辨率情况下每一个粗配准图像像素点的位置向量其中：K为当前分辨率级别，

为通过偏移公示计算得到的偏移变换，N为迭代次数；

3.3)计算两幅图像间的互信息对当前偏移变换

的梯度，当两幅图像间的互信息的增量小于预设阈值时则认为当前分辨率级数下的迭代已经收敛，将分辨率级数K递增1并进入步骤3.4)；否则返回步骤3.2)且迭代次数为N+1，进行下一次迭代，直至两幅图像间的互信息的增量不小于互信息阈值；

3.4)对前一级分辨率级数下得到的偏移变换

进行重采样，得到当前分辨率级数下的初始偏移变换，返回步骤3.2)，并进行当前分辨率级数下的迭代，直至该分辨率下两幅图像间的互信息的增量不小于互信息阈值，并且各个分辨率级数均执行迭代结束。

8.根据权利要求1或2所述的方法，其特征是，所述的插值是指：采用线性插值和中心差分方法对粗配准图像进行插值。

9.根据权利要求1或2所述的方法，其特征是，所述的精配准图像通过相关性函数进行校核：对插值后的图像以相关系数为相似性测度进行校核，设插值后图像A(x，y)和参考图像B(x，y)，相关性函数为： $C_{A,B}\left(T\right)=\frac{{\mathrm{\&Sigma;}}_{x_1,y_1}^{x_n,y_n}A\left(T(x,y)\right)\&CenterDot;B(x,y)}{\sqrt{{\mathrm{\&Sigma;}}_{x_1,y_1}^{x_n,y_n}{A}^2\left(T(x,y)\right)}\&CenterDot;\sqrt{{\mathrm{\&Sigma;}}_{x_1,y_1}^{x_n,y_n}{B}^2(x,y)}},$ 其中：A(x，y)和B(x，y)代表图像在(x，y)坐标处的灰度值，T代表图像A和图像B坐标变换，当相关系数C_A，B(T)大于等于相关系数阈值则校核认为达到最佳配准。

10.根据权利要求9所述的方法，其特征是，所述的相关系数阈值为0.985。

Images