基于广义模糊梯度矢量流场的医学序列图像运动估计方法

摘要

本发明公开了一种基于广义模糊梯度矢量流场的医学序列图像运动估计方法，包括以下步骤：1.获取一个序列图像；2.获得双步跟踪模型的广义模糊梯度矢量流场、广义模糊梯度矢量流场的局部相关性和光流矢量场；3.用手工在勾勒出第一帧图象的感兴趣区边缘轮廓；4.在三种外力场的作用下，用双步跟踪模型逐帧跟踪勾勒的感兴趣区边缘轮廓；5.结合上述轮廓线跟踪结果，用最大后验估计对轮廓线上的每一点进行优化估计与跟踪，由此得到点的最佳运动轨迹。本发明能从根本解决上述梯度矢量流场(GVF)外力场所遇到的问题，完成动态轮廓线的鲁棒跟踪并进一步实现逐点的估计与优化。

说明书

技术领域

本发明涉及图象处理方法，尤其是涉及一种反映人体心肺脏器和血管收缩、舒张运动的医学序列图像运动的估计方法。

背景技术

在医学图像后处理及图像引导的计算机辅助外科手术治疗(IGS)等领域中，基于医学图像序列的生物软组织形变和运动估计是一个重要的研究内容。

计算机视觉领域的运动估计和跟踪问题的全面研究始于上世纪八十年代初。经过不断的深入发展，已涌现出许多涉及非刚体内部的点、曲线段、轮廓以及表面等等的运动状态估计与跟踪技术。该领域中，描述序列图像中感兴趣区域、边缘、或轮廓线的时空运动状态具有重要的现实意义，各种用来计算优化动态轮廓线的方法主要有：有限差分法、有限元法、动态规划法、模拟退火法等等。用动态轮廓线模型(ACM或Snake)方法进行感兴趣区域的分割是一种较成熟、有效并且应用广泛的方法，同时，描述序列图像中动态轮廓线的运动和变化情况在医学图像后处理与计算机辅助诊断中具有更重要的研究应用价值。动态轮廓线模型这一重要分析手段的基本性质是：能对单幅图像感兴趣的区域(ROI)做局部的分割和搜寻；当施加合理的内力、外力平衡条件后，它会稳定在目标区域，形成闭合链码。如果直接将传统的动态轮廓线模型方法用于非刚体运动跟踪，其鲁棒性面临两大挑战：其一，如果动态轮廓线缺少足够的动态范围(比如离真实边界太远)，会逼近伪目标；其二，动态轮廓线对序列图像中的伪边、实边缺少足够的辨别能力。总而言之，如何产生新的外力条件、如何利用动态轮廓线模型造鲁棒的运动跟踪似然模型是该领域长期未能彻底解决的问题。本发明主要参考如下四篇文献，并以此为基础进行了改进与创新：

【1】江浩等“基于一种轮廓线估计惯性Snake模型的一般视频跟踪方法”.纽约图象处理国际论坛论文集，2002年9月22-25，页码：1301-1305。

【2】意法拉.米奇等“超声心动扫描序列图象的分割与跟踪技术：光流估计引导的动态轮廓线模型”.IEEE，医学成像，1998年，17卷(2期)，页码：127-136。

【3】许辰阳等“梯度矢量流形变模型”.霍普金斯大学学术出版社2000年9月出版的医学成像手册。

【4】陈武凡，鲁贤庆，“彩色图象边缘检测的新算法广义模糊算子法”。中国科学(A辑)，1995年，25卷(2期)，页码：219～224。

依据图像本身的性质提高轮廓线的动态范围并且优化构造曲线的外部作用力是轮廓跟踪的关键，也是完整实现非刚体运动跟踪的第一步。文献【1】公开将光流场和帧间局部相关性作为动态轮廓线的两种外力，并成功地解决了非医学视频图像的运动跟踪问题。然而，这种帧间局部相关性应用于心脏非刚体运动估计时，往往以计算量大且相关性较弱而告失败。文献【2】中公开借助光流场用有限元法成功地在超声心动序列图像中完成了大运动区域跟踪，且避免了伪边界的干扰。文献【3】中公开梯度矢量流场(GVF)作为新的外力条件在单幅图像中约束动态轮廓线，这样一来，不仅初始动态轮廓线的选取可以有更大的动态范围，而且能够逼近纯梯度场所不能达到的边缘凹陷区域，然而，利用梯度矢量流外力场分析心脏的单帧感兴趣区域时，经常遇到图像中强边缘吸引并消弱了弱边缘的梯度矢量流场，实际上感兴趣区的边界经常在弱边缘处，从而产生了较大的跟踪误差。文献【4】所提出的广义模糊理论及其边缘提取方法为本文的动态轮廓线跟踪提供了很好的边缘选择依据和鲁棒性标准。

发明内容

本发明的目的在于提供一种基于广义模糊梯度矢量流场的医学序列图像运动估计方法，能从根本解决上述梯度矢量流场(GVF)外力场所遇到的问题，完成动态轮廓线的鲁棒跟踪并进一步实现逐点的估计与优化。

为实现上述目的，本发明包括以下步骤：

1、获取一个心动周期下的不低于20帧的连续的心脏MR和CT序列图像，并将观察部位按照适当尺寸进行截取放大；

2、获得双步跟踪模型的三种外力场：第一种外力场为反映了单帧图像内部各点的空间相关性的广义模糊梯度矢量流场；第二种外力场为帧间动态轮廓线上各点周围的广义模糊梯度矢量流场的局部相关性；第三种外力场为反映了图像帧之间各点的运动相关性的光流矢量场；

3、针对第1步获得的截取放大图象，用手工勾勒出第一帧图象的感兴趣区轮廓；

4、在第2步所得外力场的作用下，用双步跟踪模型逐帧跟踪第3步得到的感兴趣区轮廓：

5、结合上述轮廓线跟踪结果，用最大后验估计对轮廓线上的每一点进行优化估计与跟踪，由此得到点的最佳运动轨迹。

本发明步骤2中获取广义模糊梯度矢量流场的具体步骤为：

a、针对步骤1中截取后的图像，逐帧获取它的广义模糊边缘图并获取其梯度；

b、分别利用平滑项自适应系数和数据项自适应系数替换经典梯度矢量流场扩散方程中的常数系数构造广义模糊梯度矢量流场扩散方程；

c、利用上述构造的广义模糊梯度矢量流场扩散方程逐帧计算图象的广义模糊梯度矢量流场

本发明步骤4中用双步跟踪模型逐帧跟踪感兴趣区轮廓线，具体步骤为：

a、第1帧初始链码的静态逼近与获取第1帧轮廓线的广义模糊梯度矢量流局部相关性外力条件：

从第1帧图像中选定感兴趣区域初始轮廓(用手工绘出)，此时便立即得到其初始链码，该初始链码在第1帧图象的广义模糊梯度矢量流场的作用下，通过双步跟踪模型的静态跟踪算子的计算，会逼近感兴趣区的真实轮廓，并形成链码，并称之为收敛态，此时第1帧轮廓跟踪完毕；

然后，利用已获得的各帧广义模糊梯度矢量流场数据，并依据已有的局部相关性算法，获得上述第1帧收敛态链码上各点的广义模糊梯度矢量流场局部相关性外力数据；

b、第2帧的动态估计：

当已经跟踪完第1帧，并且产生收敛态，便要进行第2帧的动态跟踪；此时首先进行状态赋值，即把上述步骤a中得到的第1帧的收敛态链码作为第2帧的预估计态链码；然后，在第1帧光流场和上述步骤a中得到的广义模糊梯度矢量流场局部相关性两种外力的作用下，通过双步跟踪模型的动态跟踪算子的计算，产生第2帧的估计态链码；

c、第2帧的静态逼近与获取第2帧轮廓线的广义模糊梯度矢量流局部相关性外力条件：

与第一步情况相同：首先状态赋值，即把上述步骤b中得到的第2帧的估计态链码作为第2帧的预收敛态链码；然后，在第2帧的广义模糊梯度矢量流的作用下，通过双步跟踪模型的静态跟踪算子的计算，产生第2帧的收敛态链码，此时完成了从第1帧到第2帧的跟踪过程；

然后，利用已获得的各帧广义模糊梯度矢量流数据，并依据已有的局部相关性算法，获得第2帧的收敛态链码上各点的广义模糊梯度矢量流局部相关性外力数据；

d、重复步骤b、c中的处理过程，直到最后一帧图象跟踪完毕。

本发明所述步骤5中的逐点优化估计与跟踪的具体步骤：

a、给出起始点的初始分布；

b、对初始分布的每一点，产生若干个试探点，在其中找出其与下帧动态轮廓线的最近点；

c、针对上述步骤b中的所有试探点以空间一致性和时间连续性的要求构造先验函数；

d、针对上述步骤b中的所有试探点和其对应下帧动态轮廓线上的最近点，构造似然约束条件，并获得似然概率；

e、利用上述构造的先验约束条件和似然约束条件，对上述具有马尔可夫随机特性的试探点进行最大后验估计、并逐帧获得最大后验概率，可得到特征点的最佳运动轨迹，从而完成感兴趣区的动态跟踪。

对本发明的比较试验如下：1、请心外科专家对两类心脏图像逐帧勾勒出心室和心房形变运动轮廓(如图3)；2、选择每类图像的第一帧所勾勒出的轮廓作为初始链码，用本发明进行逐帧处理；3、分别用梯度矢量流场与广义模糊梯度矢量流场两种外力作用于双步跟踪模型进行跟踪，将两类跟踪结果与手工勾勒出的结果相对比既有直观的差别(如图4)，更有量化均方误差(见表1)，由表可见广义模糊梯度矢量流(GFGVF)场外力条件下的跟踪精度明显好于梯度矢量流(GFGVF)场条件。

                [表1]跟踪结果误差对比

因此采用广义模糊梯度矢量流场，使得梯度流场得以优化，图像平坦处的梯度流数据得以更好的平滑、图像边缘处的梯度数据得以更好地恢复，从而解决了全局和局部适应性的矛盾；避免了用经典梯度矢量流场处理心脏图像序列时，多帧图像出现汇流外溢而导致动态轮廓线出现异常的形变结果，提高了鲁棒性。

附图说明

图1为本发明的流程框图；

图2为单个周期MR图像序列下的心脏左心室内壁形变跟踪过程描述。这种基于广义模糊梯度矢量流场的双步跟踪模型可对序列图像的感兴趣区的轮廓进行鲁棒的时空跟踪。

图3为针对CT六帧心脏序列图像的左心房，心外科医生手工描绘的感兴趣的轮廓边缘；

图4为针对图3中CT六帧心脏序列图像的左心房，分别将梯度矢量流场(上行)与广义模糊梯度矢量流场(下行)两种外力作用下的双步跟踪模型模型轮廓跟踪结果与图3心外科医生手工描绘边缘进行对比。比较结果不仅显示出双步跟踪模型跟踪算法的稳定性，同时也说明广义模糊梯度矢量流场具有更优的性质。

具体实施方式

1、获取一个心动周期下的连续的心脏MR和CT序列图像25帧，并利用现有的插值算法将观察部位按照适当尺寸进行截取放大，该方法有利于增强图像感兴趣区的细节分辨率、并提高跟踪的质量；

2、利用已有的光流计算方法逐帧计算步骤1中截取放大后图像的光流场；

3、针对截取后的图像，利用现有方法逐帧计算它的广义模糊边缘图；

4、获取广义模糊梯度矢量流场，具体步骤如下：

a、获取步骤3中所得的广义模糊图象的广义模糊边缘数据I^e并计算其梯度I^e；

b、构造广义模糊梯度矢量流场扩散方程：分别利用平滑项自适应系数和数据项自适应系数ηexp(-(|μ′_I|/σ)²)和ρ(1-g(·))|I^e|²替换经典梯度矢量流场扩散方程中的常数系数η和|I^e|²；构造出的广义模糊梯度矢量流场扩散方

程为：

          U_t＝g(|μ′^I|)²U-ρ(1-g(μ′_I))| I^e|²(U-I^e)

c、利用如上广义模糊梯度矢量流方程逐帧计算图象的广义模糊梯度矢量流场；

5、针对第1步获得的截取放大图象，手工勾勒出第一帧图象的感兴趣区轮廓；

6、用双步跟踪模型逐帧跟踪步骤5中的动态轮廓线，双步跟踪模型的理论为：双步跟踪模型由静态算子和动态算子构成，它们分别用来解决运动的静态逼近和动态估计的问题。其中动态算子为经典梯度矢量流模型的改进，在经典梯度矢量流模型的算法的结构形式中改变了经典梯度矢量流模型的受力条件，即增加了光流矢量场和动态轮廓线上各点的广义模糊梯度矢量流的局部相关性两种外力；其静态算子为经典梯度矢量流模型的算法的结构形式，但采用了广义模糊梯度矢量流外力场而不是其原本的梯度矢量流外力场。逐帧跟踪动态轮廓线具体步骤如下：

a、第1帧初始链码的静态逼近与获取第1帧轮廓线的广义模糊梯度矢量流场局部相关性外力条件：

从第1帧图像中选定感兴趣区域初始轮廓(用手工绘出)，此时便立即得到其初始链码，该初始链码在第1帧图象的广义模糊梯度矢量流场的作用下，通过双步跟踪模型的静态跟踪算子的计算，会逼近感兴趣区的真实轮廓，并形成链码，并称之为收敛态，此时第1帧轮廓跟踪完毕；

然后，利用已获得的各帧广义模糊梯度矢量流场数据，并依据已有的局部相关性算法，获得上述第1帧收敛态链码上各点的广义模糊梯度矢量流场局部相关性外力数据；

b、第2帧的动态估计：

当已经跟踪完第1帧，并且产生收敛态，便要进行第2帧的动态跟踪；此时首先进行状态赋值，即把上述步骤a中得到的第1帧的收敛态链码作为第2帧的预估计态链码；然后，在第1帧光流场和上述步骤a中得到的广义模糊梯度矢量流场局部相关性两种外力的作用下，通过双步跟踪模型的动态跟踪算子的计算，产生第2帧的估计态链码；

c、第2帧的静态逼近与获取第2帧轮廓线的广义模糊梯度矢量流场局部相关性外力条件：

与第一步情况相同：首先状态赋值，即即把上述步骤b中得到的第2帧的估计态链码作为第2帧的预收敛态链码；然后，在第2帧的广义模糊梯度矢量流场的作用下，通过双步跟踪模型的静态跟踪算子的计算，产生第2帧的收敛态链码，此时完成了从第1帧到第2帧的跟踪过程；

然后，利用已获得的各帧广义模糊梯度矢量流场数据，并依据已有的局部相关性算法，获得第2帧的收敛态链码上各点的广义模糊梯度矢量流场局部相关性外力数据；

d、重复步骤b、c中的处理过程，直到最后一帧图象跟踪完毕；

7、结合上述轮廓线跟踪结果，用最大后验估计对轮廓线上的每一点进行优化估计与跟踪，具体步骤：

a、给出起始点的初始分布；

b、对初始分布的每一点，产生64个试探点，在其中找出其与下帧动态轮廓线的最近点；试探点越多越好，但太多会增加计算的复杂度；

c、针对上述步骤b中的所有试探点以空间一致性和时间连续性的要求构造先验函数；上述空间一致性和时间连续性的理论为：被视为运动整体的刚性和非刚性物体中的某一质点(或微粒)具有和其邻点一致或相近的运动状态，该运动状态随时间进行连续变化。空间一致性和时间连续性可用贝叶斯方法表示为如下概率形式：如果n代表具有时间意义的帧数、i代表具有空间意义的特征点序数，则：

$$ >>P>>(>>X>n>>)>>=>P>>(sup>>x>n>isup>>|sup>>x>>n>->1>>isup>>)>>P>>(sup>>x>n>>i>+>1>sup>>|sup>>x>n>isup>>)>>>$$

上述第一项概率是时间连续性条件、第二项是空间一致性条件，以概率形式始终保持该两项乘积最大便是空间一致性和时间连续性的基本要求；

d、针对上述步骤b中的所有试探点和其对应下帧动态轮廓线上的最近点，构造似然约束条件，并获得似然概率；

e、利用上述构造的先验约束条件和似然约束条件，对上述具有马尔可夫随机特性的试探点进行最大后验估计、并逐帧获得最大后验概率，可得到每一点的最佳运动轨迹，从而完成感兴趣区的动态跟踪。