基于运动预测和三维约束的序列图像分割方法

摘要

本发明公开了一种基于运动预测和偏三维约束并结合主动轮廓模型的图像序列分割方法，属于图像处理技术领域。以主动轮廓模型上的控制点作为预测模型形变的着手点，把物体变化的平滑性和一致性，表示为控制点位置变化的一致性，以控制点的变化趋势，预测在序列中下一帧图像中相应控制点的位置，然后以预测位置作为匹配的初始搜索位置采用MPEG－4所使用的是菱形搜索算法。本发明引入了偏三维约束的概念，作为主动轮廓模型的能量之一。这样又尽量避免了错误匹配造成的影响。达到了在简单、计算量小的条件下的较好的预测效果。

说明书

技术领域

本发明属于图像处理技术领域，具体涉及一种基于运动预测和偏三维约束，并结合主动轮廓模型的图像序列分割方法。

背景技术

三维图像分割是图像处理和计算机视觉中的一个重要研究领域，广泛应用于计算机辅助三维重建、运动物体跟踪以及医学中细胞运动和形变的分析等应用领域。三维图像或图像序列与一般图像的不同之处在于：序列图像是针对同一目标，切片是由间隔相同的距离或时间连续采样而得到的，因此对三维图像的分割既与一般的图像分割有相同之处，又有特殊之处。

传统的主动轮廓模型都要求模型初始化在目标轮廓的附近，否则将可能收敛到错误的结果。而在三维图像分割中对每一幅切片图进行手工初始化是一种费时而繁琐的工作，需要找到一种自动化的三维分割方法。

现有的基于主动轮廓模型的分割方法主要有：

1、预测法，预测法简单易行，利用物体的轮廓一般都是平滑变化这一常识，用前面数个轮廓的位置，对初始位置的设置进行改进，在自动分割中，如何更精确的自动设定主动轮廓的初始位置非常重要，初始位置距离实际边缘越接近，收敛到正确位置的可能性就越大。这种方法有一定效果，但没有利用图像本身的信息，有进一步改进的空间。

2、光流法，光流法利用光流场的方法迭代获得图像中物体运动的趋势，这种趋势以一个新的外力的形式影响收敛的结果，有较好的效果。但缺点是方法复杂，运算量比较大。而且没有改善初始化位置，迭代次数也会比较多，会增加收敛的时间。

3、三维形变模型法，常规的三维形变模型法是把物体整体建立三维轮廓模型，虽然利用了相邻图像相关性信息但是模型复杂，计算量大，初始位置设置困难。如果能在二维轮廓模型中引入三维约束，则可以在二维模型的优点中加入三维模型的三维平滑的特点，更符合物体的实际特性。

发明内容

技术问题：本发明基于三维图像序列的相似性，采用在视频压缩中常见的运动估计方法，并与偏三维约束相结合，提出了一种新的适用于三维图像分割的方法。基于三维图像具有三维平滑、连续变化等特性，本发明提出了一种新的图像序列分割方法。该方法在计算过程中只需手动设置第一幅图中主动形变模型的关键点位置。在其他图像中，首先采用预测和视频运动估计中常用的块匹配法优化，在MPEG-4和H.26x中采用的都是这种块匹配法，然后采用偏三维约束和梯度矢量流(GVF)，从初始化位置开始在图像中进一步迭代收敛得到每幅图像中的最终轮廓位置，在一般的二维平滑基础上，达到三维平滑，三维分割的作用。本发明方法可以显著提高分割的准确性以及速度。

技术方案：本发明公开了一种基于运动预测和三维约束的序列图像分割方法，该方法包括如下步骤：

1)对数字图像初始化：

11)在第一帧数字图像轮廓线转向的地方，交互式设置手动控制点，相邻的两个手动控制点之间均等生成控制点的初始位置；

12)迭代收敛，收敛到第一幅图像中物体轮廓的边沿，其过程如下：

121)采取加入控制点或合并控制点的方式来调节相邻两个控制点的间距，

122)由蛇的控制点的初始位置出发，以梯度矢量流场替代原图像计算梯度值并作为外力，

123)计算第一层内能公式为E_int＝(x_i+1，j+x_i-1，j-2x_i，j)，x_i，j表示第j幅图像中的第i个控制点的位置，E_int为内部能量，

124)计算合力公式E_snake＝E_int+E_image+E_constrain，计算出的合力为零来作为图像收敛的条件，其中E_snake为能量函数，E_constrain为约束能量来限制控制点的运动，E_int为内部能量，E_image为外部能量，

125)根据公式$x_{i,j}=x_{i,j}+E_{\mathrm{snake}}^{i,j}$计算出来的合力移动控制点，x_i，j表示第j幅图像中的第i个控制点的位置，E_snake^i，j表示第j幅图像中的第i个控制点的位置的能量函数，

126)对图像收敛的处理过程：

利用公式：

avd＝sum(abs(E_snake))/length(E_snake)来判断图像收敛条件，sum(abs(E_snake))表示为所有控制点合力的绝对值之和，length(E_snake)为控制点的个数，avd为合力的平均值，

如果合力的平均值小于0.08时再收敛15次，目的是在合力较小的情况下进一步接近轮廓，

当avd＜0.08后且又迭代15次时，迭代结束，得到最终收敛结果，

否则，回到步骤121)，迭代继续进行；

13)对第一帧数字图像不进行控制点位置预测计算，以第一帧数字图像主动轮廓模型的收敛位置作为控制点的预测位置，

按合力公式E_snake＝E_int+E_image+E_constrain计算出的合力为零来作为图像收敛的条件，其中E_snake为能量函数，E_constrain为约束能量来限制控制点的运动，E_int为内部能量，E_image为外部能量；

14)块匹配初始搜索点的选择，以预测的下一帧的控制点的位置作为初始搜索点，和直接采用参考图像对应位置相比，2个图像块比较，相应的像素点两两求差的绝对值，它们的平均值是最小绝对差匹配准则所求的值，公式中的M，N是图像块的长和宽，f_k(m，n)是指在k帧数字图像中(m，n)处像素点的像素值，

${\mathrm{MAD}}_{(\mathrm{\Delta p},\mathrm{\Delta q})}=(1/\mathrm{MN})\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}|f_k(m,n)-f_{k-1}(m+\mathrm{\Delta p},n+\mathrm{\Delta q})|$

其中(Δp，Δq)为位移矢量，f_k为当前帧图像的灰度值，f_k-1为上一帧图像的灰度值，图像块大小为M*N，MAD_(Δp，Δq)为最小值是最优匹配点，

搜索策略为菱形搜索算法，采用9个检测点的大模板和5个检测点的小模板，对匹配得到的平均绝对差设定门限值，当最终搜索到最优匹配点小于门限值时，把匹配到的点作为控制点的初始位置；反之，直接把预测位置作为控制点的初始位置；

2)循环：

21)设置数字图像帧层数记数器，用于判断是否超过最底数字图像帧层，若数字图像帧层数没有超过最底数字图像帧层则执行步骤22)，否则结束循环；

22)以前两层收敛的控制点的坐标预测下一层的控制点的位置

预测下一帧数字图像的控制点位置，x_i，j表示第j帧数字图像中的第i个控制点的位置，控制点位置的变化用V_i，j表示为：V_i，j＝x_i，j-x_i，j-1，得到预测的j+1帧数字图像中的相应控制点的位置：${\tilde{x}}_{i,j+1}=x_{i,j}+k_1\times V_{i,j},$k₁为预测系数，当k₁为0时，直接用上一帧数字图像的分割结果作为初始值；

23)以步骤22)中所预测控制点的位置为起点作为初始搜索点，和直接采用参考图像对应位置相比，2个图像块比较，相应的像素点两两求差的绝对值，它们的平均值是最小绝对差匹配准则所求的值，公式中的M，N是图像块的长和宽，f_k(m，n)是指在k帧数字图像中(m，n)处像素点的像素值，

${\mathrm{MAD}}_{(\mathrm{\Delta p},\mathrm{\Delta q})}=(1/\mathrm{MN})\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}|f_k(m,n)-f_{k-1}(m+\mathrm{\Delta p},n+\mathrm{\Delta q})|$

其中(Δp，Δq)为位移矢量，f_k为当前帧图像的灰度值，f_k-1为上一帧图像的灰度值，图像块大小为M*N，MAD_(Δp，Δq)为最小值是最优匹配点，

搜索策略为菱形搜索算法，采用9个检测点的大模板和5个检测点的小模板，对匹配得到的平均绝对差设定门限值，当最终搜索到最优匹配点小于门限值时，把匹配到的点作为控制点的初始位置；反之，直接把预测位置作为控制点的初始位置；

24)偏三维约束的迭代收敛，过程如下：

241)采取加入控制点或合并控制点的方式来调节相邻两个控制点的间距，

242)由步骤22)中所预测控制点的位置出发，以梯度矢量流场替代原图像计算梯度值并作为外力，

243)计算偏三维约束的内力公式为

E_int′＝k(x_i+1，j+x_i-1，j-2x_i，j)+(1-k)(x_i，j+1+k₁(x_i，j+1-x_i，j+2)-x_i，j)，

x_i，j表示第j幅图像中的第i个控制点的位置，E_int为偏三维约束内部能量，k为预测系数，

244)合力公式为E_snake＝E_int+E_image+E_constrain，计算出的合力为零来作为图像收敛的条件，其中E_snake为能量函数，E_constrain为约束能量来限制控制点的运动，E_int内部能量和E_image外部能量；

245)合力移动控制点，$x_{i,j}=x_{i,j}+E_{\mathrm{snake}}^{i,j},$x_i，j表示第j幅图像中的第i个控制点的位置，E_snake^i，j表示第j幅图像中的第i个控制点的位置的能量函数，

246)判断是否到达收敛，

利用公式：

avd＝sum(abs(E_snake))/length(E_snake)作为判断图像收敛条件，sum(abs(E_snake))表示为所有控制点合力的绝对值之和，length(E_snake)为控制点的个数，avd为合力的平均值，

如果合力的平均值小于0.08时再收敛15次，目的是在合力较小的情况下进一步接近轮廓，

当avd＜0.08后且又迭代15次时，迭代结束，得到最终收敛结果，

否则，回到步骤241)，迭代继续进行；

25)数字图像帧层数记数器进行帧层数累加；

26)判断步骤25)中的数字图像帧层数是否超过最底数字图像帧层，若数字图像帧层数没有超过最底数字图像帧层则执行步骤22)，否则结束循环。

上述基于运动预测和三维约束的序列图像分割方法，该方法步骤14)和步骤23)中的门限值为10到20之间，步骤121)和步骤241)中的两个控制点的间距保持在3个像素到7个像素之间，步骤243)中的k1的预测系数取值为0到1之间。

本发明首先采用预测与匹配相结合初始位置设置方法。以主动轮廓模型上的控制点作为预测模型形变的着手点，避免了有的预测法中关于质心的概念。对质心的应用中要求轮廓为比较饱满的形状，否则若质心太过接近于轮廓边缘，会使判断的法线方向有较大的偏差，影响有效性。我们把物体变化的平滑性和一致性，表示为控制点位置变化的一致性。

由于我们是以控制点构成的网格表示物体的轮廓，物体表现出来的平滑性和一致性要体现在控制点网格上以控制点的变化趋势，预测在序列中下一帧图像中相应控制点的位置，然后以预测位置作为匹配的初始搜索位置采用MPEG-4所使用的是菱形搜索算法。对于匹配结果，由于控制点是位于物体的边缘部分，这个区域图像灰度值变化大，有利于匹配到正确的结果。并进一步设定一个恰当门限，图像变化以及噪声剧烈的，门限就要高一点，变化及噪声平稳的就小一点。当最终搜索到最优匹配点的平均绝对差(MAD)小于门限时，就认为结果为正确匹配点，把匹配到的点作为控制点的初始位置；反之则认为结果为错误的匹配点，忽略匹配结果，直接把预测位置作为控制点的初始位置。这样又尽量避免了错误匹配造成的影响。达到了在简单，计算量小的条件下的较好的效果。

在传统的主动轮廓模型基础上，本发明引入了偏三维约束的概念，作为主动轮廓模型的能量之一。对于本发明的三维自动分割方法，是从第一层图像开始，逐层向下。在任何时候，我们只知道上一层图像中的控制点的位置，而不知道下一层图像中的控制点位置，所以并不能直接采用4点约束式的三维约束模型。这个种情况下，本发明基于上层控制点和预测，提出了偏三维约束。

对于横向内能，目的是使轮廓在水平方向上平滑，由同层相邻的两个控制点决定，和三维模型没有区别。而纵向内能，由于本文序列分割的特性，是由二维模拟三维，不能轻易得出，因为纵向相邻的2个控制点只知道一个。对此，本发明采用预测的方法来进行模拟。由于物体的平滑性，物体在每层图像上投影的变化具有连续性，即由上两层图像控制点预测得到的位置与上下两层控制点的中点的位置在大多数情况下近似，并且同样能反映物体平滑的特性。

通过手动设置第一幅图中主动形变模型的关键点位置后，在其他图像中，采用预测和块匹配法相结合的方法自动地设置后续图像的主动轮廓模型初始化位置，并采用偏三维约束和梯度矢量流相结合的方法，从初始位置开始，自动化地进一步迭代收敛得到每幅图像中的最终轮廓位置。

采用预测与匹配相结合初始位置设置方法。以主动轮廓模型上的控制点作为预测模型形变的着手点，把物体变化的平滑性和一致性，表示为控制点位置变化的一致性，以控制点的变化趋势，预测在序列中下一帧图像中相应控制点的位置，然后以预测位置作为匹配的初始搜索位置采用MPEG-4所使用的是菱形搜索算法。对于匹配结果，由于控制点是位于物体的边缘部分，这个区域图像灰度值变化大，有利于匹配到正确的结果。并进一步设定一个恰当门限。当最终搜索到最优匹配点的MAD小于门限时，就认为结果为正确匹配点，把匹配到的点作为控制点的初始位置；反之则认为结果为错误的匹配点，忽略匹配结果，直接把预测位置作为控制点的初始位置。这样又尽量避免了错误匹配造成的影响。达到了在简单、计算量小的条件下的较好的预测效果。

在传统的主动轮廓模型基础上，本发明引入了偏三维约束的概念，作为主动轮廓模型的能量之一。对于本发明的三维自动分割方法，是从第一层图像开始，逐层向下。对于横向内能，由同层相邻的两个控制点决定，和三维模型没有区别。而纵向内能则以上层图像中的控制点的位置采用偏三维约束，以预测的方法来进行模拟。

有益效果：本发明中联合匹配预测用于序列图像分割，优化了初始主动轮廓的设置，使得初始轮廓更接近与实际边缘，有利于得到准确的分割结果。在二维轮廓模型中引入三维约束，则可以在二维模型的优点中加入三维模型的三维平滑的特点，更符合物体的实际特性，并可以同时实现重建的作用，以低计算成本达到较好的效果。该方法在三维序列图像分割中非常有效，可以显著提高分割的准确性以及速度。

附图说明

图1为本发明总体流程图。

图2为偏三维约束迭代收敛过程的流程图。

图3为预测法对人体肝脏CT图像序列进行分割的成像结果。

图4为匹配法对人体肝脏CT图像序列进行分割的成像结果。

图5为光流法对人体肝脏CT图像序列进行分割的成像结果。

图6为本发明方法对人体肝脏CT图像序列进行分割的成像结果。

具体实施方式

下面是本发明的具体实施例来进一步描述：

本发明的运行环境是WinXP操作系统；闪龙2500+CPU；算法采用仿真计算平台Matlab 7.0版本进行编制。该方法总体过程如图1所示包括如下步骤：

1)对数字图像初始化：

11)在第一帧数字图像轮廓线转向的地方，交互式设置手动控制点，相邻的两个手动控制点之间均等生成控制点的初始位置；

12)迭代收敛，收敛到第一幅图像中物体轮廓的边沿，包括如下步骤：

121)采取加入控制点或合并控制点的方式来调节相邻两个控制点的间距保持在5个像素，

122)由蛇的控制点的初始位置出发，以梯度矢量流场替代原图像计算梯度值并作为外力，

123)计算第一层内能公式为E_int＝(x_i+1，j+x_i-1，j-2x_i，j)，x_i，j表示第j幅图像中的第i个控制点的位置，E_int为内部能量，

124)计算合力公式E_snake＝E_int+E_image+E_constrain，计算出的合力为零来作为图像收敛的条件，其中E_snake为能量函数，E_constrain为约束能量来限制控制点的运动，E_int为内部能量，E_image为外部能量，

125)根据公式$x_{i,j}=x_{i,j}+E_{\mathrm{snake}}^{i,j}$计算出来的合力移动控制点，x_i，j表示第j幅图像中的第i个控制点的位置，E_snake^i，j表示第j幅图像中的第i个控制点的位置的能量函数，

126)对图像收敛的处理过程：

利用公式：

avd＝sum(abs(E_snake))/length(E_snake)来判断图像收敛条件，sum(abs(E_snake))表示为所有控制点合力的绝对值之和，length(E_snake)为控制点的个数，avd为合力的平均值，

如果合力的平均值小于0.08时再收敛15次，目的是在合力较小的情况下进一步接近轮廓，

当avd＜0.08后且又迭代15次时，迭代结束，得到最终收敛结果，

否则，回到步骤121)，迭代继续进行；

13)对第一帧数字图像不进行控制点位置预测计算，以第一帧数字图像主动轮廓模型的收敛位置作为控制点的预测位置，

按合力公式E_snake＝E_int+E_image+E_constrain计算出的合力为零来作为图像收敛的条件，其中E_snake为能量函数，E_constrain为约束能量来限制控制点的运动，E_int为内部能量，E_image为外部能量；

14)块匹配初始搜索点的选择，以预测的下一帧的控制点的位置作为初始搜索点，和直接采用参考图像对应位置相比，2个图像块比较，相应的像素点两两求差的绝对值，它们的平均值是最小绝对差匹配准则所求的值，公式中的M，N是图像块的长和宽，f_k(m，n)是指在k帧数字图像中(m，n)处像素点的像素值，

${\mathrm{MAD}}_{(\mathrm{\Delta p},\mathrm{\Delta q})}=(1/\mathrm{MN})\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}|f_k(m,n)-f_{k-1}(m+\mathrm{\Delta p},n+\mathrm{\Delta q})|$

其中(Δp，Δq)为位移矢量，f_k为当前帧图像的灰度值，f_k-1为上一帧图像的灰度值，图像块大小为M*N，MAD_(Δp，Δq)为最小值是最优匹配点，

搜索策略为菱形搜索算法，采用9个检测点的大模板和5个检测点的小模板，对匹配得到的平均绝对差设定门限值为15，当最终搜索到最优匹配点小于门限值时，把匹配到的点作为控制点的初始位置；反之，直接把预测位置作为控制点的初始位置；

2)循环：

21)设置数字图像帧层数记数器，用于判断是否超过最底数字图像帧层，若数字图像帧层数没有超过最底数字图像帧层则执行步骤22)，否则结束循环；

22)以前两层收敛的控制点的坐标预测下一层的控制点的位置

预测下一帧数字图像的控制点位置，x_i，j表示第j帧数字图像中的第i个控制点的位置，控制点位置的变化用V_i，j表示为：V_i，j＝x_i，j-x_i，j-1，得到预测的j+1帧数字图像中的相应控制点的位置：${\tilde{x}}_{i,j+1}=x_{i,j}+k_1\times V_{i,j},$k₁为预测系数，当k₁为0时，直接用上一帧数字图像的分割结果作为初始值；

23)以步骤22)中所预测控制点的位置为起点作为初始搜索点，和直接采用参考图像对应位置相比，2个图像块比较，相应的像素点两两求差的绝对值，它们的平均值是最小绝对差匹配准则所求的值，公式中的M，N是图像块的长和宽，f_k(m，n)是指在k帧数字图像中(m，n)处像素点的像素值，

${\mathrm{MAD}}_{(\mathrm{\Delta p},\mathrm{\Delta q})}=(1/\mathrm{MN})\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}|f_k(m,n)-f_{k-1}(m+\mathrm{\Delta p},n+\mathrm{\Delta q})|$

其中(Δp，Δq)为位移矢量，f_k为当前帧图像的灰度值，f_k-1为上一帧图像的灰度值，图像块大小为M*N，MAD_(Δp，Δq)为最小值是最优匹配点，

搜索策略为菱形搜索算法，采用9个检测点的大模板和5个检测点的小模板，对匹配得到的平均绝对差设定门限值，当最终搜索到最优匹配点小于门限值时，把匹配到的点作为控制点的初始位置；反之，直接把预测位置作为控制点的初始位置；

24)偏三维约束的迭代收敛，其过程如图2所示包括如下步骤：

241)采取加入控制点或合并控制点的方式来调节相邻两个控制点的间距，

242)由步骤22)中所预测控制点的位置出发，以梯度矢量流场替代原图像计算梯度值并作为外力，

243)计算偏三维约束的内力公式为

E_int′＝k(x_i+1，j+x_i-1，j-2x_i，j)+(1-k)(x_i，j+1+k₁(x_i，j+1-x_i，j+2)-x_i，j)，

x_i，j表示第j幅图像中的第i个控制点的位置，E_int为偏三维约束内部能量，k1为预测系数取0.5，

244)合力公式为E_snake＝E_int+E_image+E_constrain，计算出的合力为零来作为图像收敛的条件，其中E_snake为能量函数，E_constrain为约束能量来限制控制点的运动，E_int内部能量和E_image外部能量；

245)合力移动控制点，$x_{i,j}=x_{i,j}+E_{\mathrm{snake}}^{i,j},$x_i，j表示第j幅图像中的第i个控制点的位置，E_snake^i，j表示第j幅图像中的第i个控制点的位置的能量函数，

246)判断是否到达收敛，

利用公式：

avd＝sum(abs(E_snake))/length(E_snake)作为判断图像收敛条件，sum(abs(E_snake))表示为所有控制点合力的绝对值之和，length(E_snake)为控制点的个数，avd为合力的平均值，

如果合力的平均值小于0.08时再收敛15次，目的是在合力较小的情况下进一步接近轮廓，

当avd＜0.08后且又迭代15次时，迭代结束，得到最终收敛结果，

否则，回到步骤241)，迭代继续进行；

25)数字图像帧层数记数器进行帧层数累加；

26)判断步骤25)中的数字图像帧层数是否超过最底数字图像帧层，若数字图像帧层数没有超过最底数字图像帧层则执行步骤22)，否则结束循环。

在第一幅序列图像上，沿着所需的图像，手动设置是以鼠标点击的方式直接在图像边缘附近点击设定若干点，标识轮廓位置，根据手动设置的点，自动等间距的生成控制点手动设置少数的接近轮廓的点，这些点主要位于轮廓线转向的地方，在手动设置的相邻的两个点之间也许距离会比较远，而控制点的间距要求比较均匀，间距较小，根据这个原则等距的均匀地加入控制点。控制点的位置作为主动轮廓模型的初始位置。

轮廓的迭代收敛，由主动轮廓的初始位置，采用改进外力的Snake模型，GVF Snake模型收敛主动轮廓。采用公式(1)计算各个控制点的能量，以合力推动控制点的移动。公式中的变量代码能量就是指下面公式中计算出来的数值，分为x轴方向和y轴方向，这计算出来的值越大，控制点的移动也就越大。这个能量分为下面公式中的3个部分，这3个能量的和作为推动控制点移动的合力。

E_snake＝E_int+E_image+E_constrain    (1)

其中外部能量E_image采用梯度矢量流场。

离散计算公式为

$\left(\begin{array}{c}{u}^{(n+1)}(p,q)=u^{\left(n\right)}(p,q)+L{▿}^2{u}^{\left(n\right)}(p,q)-[u^{\left(n\right)}(p,q)-f_x(p,q)][f_x{(p,q)}^2+f_y{(p,q)}^2]\\ v^{(n+1)}(p,q)=v^{\left(n\right)}(p,q)+L{▿}^2{v}^{\left(n\right)}(p,q)-[v^{\left(n\right)}(p,q)-f_y(p,q)]\left[f_x{(p,q)}^2{+f}_y{(p,q)}^2\right]\end{array}\right)$

(2)

计算出梯度矢量流场，其中：

²u⁽ⁿ⁾(p，q)＝u⁽ⁿ⁾(p+1，q)+u⁽ⁿ⁾(p-1，q)+u⁽ⁿ⁾(p，q+1)+u⁽ⁿ⁾(p，q-1)+4u⁽ⁿ⁾(p，q)

²v⁽ⁿ⁾(p，q)＝v⁽ⁿ⁾(p+1，q)+v⁽ⁿ⁾(p-1，q)+v⁽ⁿ⁾(p，q+1)+v⁽ⁿ⁾(p，q-1)+4v⁽ⁿ⁾(p，q)

(3)

u⁽ⁿ⁾(p，q)，v⁽ⁿ⁾(p，q)是计算GVF场时，第n次迭代的的迭代值，第一次迭代时，它的值就是图像本身的像素值。f_x(p，q)，f_y(p，q)是x，y方向上的图像梯度值。

在三维图像中，内部能量与周围4个控制点相关，受4个控制点的共同影响，周围的4个控制点决定内部能量的大小。本发明采用预测的方法来进行模拟，采用预测位置作为纵向平滑约束，提出了偏三维约束。公式表示为：

E_int′＝k(x_i+1，j+x_i-1，j-2x_i，j)+(1-k)(x_i，j+1+k₁(x_i，j+1-x_i，j+2)-x_i，j)

(4)

其中，x_i，j表示第j幅图像中的第i个控制点的位置，k控制横向约束的和纵向预测之间的权衡，k₁为预测系数(0≤k₁≤1)。由于预测的准确性低于实际的纵向内能且上层中的控制点位置固定不变，物体的轮廓是控制点的网格，每个控制点相邻的都有4个控制点，水平方向上2个作为横向约束，产生横向能量，垂直方向上2个作为纵向约束，产生纵向内能，这是传统的三维模型计算法，本文方法在横向上是和它一致的，而在纵向上使用预测的方法产生纵向能量。约束过强是指由于预测的结果还是有一定的不准确性，为了避免这种不准确性导致错误，横向约束应强于纵向约束(0.5＜k＜1)，使轮廓有更大的自由度，避免约束过强造成的误差。在实验中，k取0.6-0.7，k₁取0.5。

在传统的三维模型中，内力的公式表示为：

E_int＝(x_i+1，j+x_i-1，j-2x_i，j)+(x_i，j+1+x_i，j-1-2x_i，j)

区别在于，控制点周围的4个控制点都是已知的，使用这4个控制点约束中间控制点。

在本算法中，第一层的约束是普通的二维收敛，采用的内能也是二维的，表现为只有横向约束，没有三维的纵向约束，公式为：

E_int＝(x_i+1，j+x_i-1，j-2x_i，j)

采用下述约束能量限制控制点的运动，防止跳变的产生控制点一次变化过大可能会落入无法回归的图像区域而而产生错误

$\bar{C}=\frac{1}{m}\underset{i=1}{\overset{m}{\Sigma }}(E_{\mathrm{int}}^i+E_{\mathrm{image}}^i)$

(5)

$\sigma =\sqrt{\frac{1}{m}\underset{i=1}{\overset{m}{\Sigma }}{(E_{\mathrm{int}}^i\text{+}{E}_{\mathrm{image}}^i-\bar{C})}^2}$

(6)

E_snake＝E_constrain(E_int+E_image)

(7)

$E_{\mathrm{constrain}}^i=\min \{1,\frac{3\sigma }{E_{\mathrm{int}}^i+E_{\mathrm{image}}^i-\bar{C}}\}$

(8)

内部能量和外部能量的合力在约束能量的约束之下控制控制点的移动。

反复计算，经过多次迭代最终收敛到所需的图像边沿处。

预测下一层的控制点位置

分两种情况处理。当为第一幅图像时，就以主动轮廓模型的收敛位置作为预测位置，不再进行预测计算。

对其它的下层图像，则进行如下的预测处理。

控制点位置的变化V_i，j可以表示为：

V_i，j＝x_i，j-x_i，j-1

(9)

x_i，j表示第j幅图像中的第i个控制点的位置，可以得到预测的n+1幅图中的相应控制点的位置：

${\tilde{x}}_{i,j+1}=x_{i,j}+k_1\times V_{i,j}$

(10)

当k₁为0时，就是直接使用上一幅图的分割结果作为初始值。

匹配。利用预测的结果进一步采用匹配的方法获得更为接近与图像中物体边缘的位置。本发明采用实现简单且被MPEG-x和H.26x等广泛应用的基于块匹配的运动估计算法块匹配是指2个相同大小的图像进行比较，差异越小匹配度越高。匹配点是指块匹配图像块的中心点，这个点可以匹配块的位置，块匹配成功实际上相应的匹配点匹配成功，块匹配技术主要包括以下3点：

第一，初始搜索点的选择。本发明中采用预测的下一帧的控制点的位置作为初始搜索点作为初始搜索点，和直接采用上一帧的控制点收敛位置作为初始搜索点相比，可以有效减少匹配的次数和时间。

第二，块匹配准则。本发明方法中采用常用的最小绝对差匹配准则2个图像块比较，相应的像素点两两求差的绝对值，它们的平均值就是最小绝对差匹配准则所求的值。公式如下：

${\mathrm{MAD}}_{(\mathrm{\Delta p},\mathrm{\Delta q})}=(1/\mathrm{MN})\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}|f_k(m,n)-f_{k-1}(m+\mathrm{\Delta p},n+\mathrm{\Delta q})|$

(11)

M，N是图像块的长和宽，f_k(m，n)是指在k帧中(m，n)处像素点的像素值，其中(Δp，Δq)为位移矢量；MAD_(Δp，Δq)的最小值是最优匹配点。

第三，搜索策略。使用的是菱形搜索算法(Diamond Search，DS)。采用了两种模板，即9个检测点的大模板LDSP(Large Diamond Search Pattern)和5个检测点的小模板SDSP(Small Diamond Search Pattern)。首先，搜索时先采用大模板对搜索中心点周围的9个点进行匹配运算，就是把以这9个点为中心的9个图像块分别和上一帧相应控制点为中心的图像块计算MAD_(Δp，Δq)，当最小块误差点即为中心点时，则自动转换为小模板搜索策略，对中心点周围其他5个点进行匹配运算，5点中的最优匹配点即为所求点。如果在采用大模板搜索时，最优点不在中心点时，那么就以该最优点为中心点重新采用大模板进行搜索，直到最优匹配点为中心点时，再转换为小模板搜索，最终找到最优匹配点。

对于匹配结果，由于控制点是位于物体的边缘部分，这个区域图像灰度值变化大，有利于匹配到正确的结果。并进一步对匹配得到的平均绝对差设定一个门限。当最终搜索到最优匹配点的MAD小于门限时，就认为结果为正确匹配点，把匹配到的点作为控制点的初始位置；反之则认为结果为错误的匹配点，忽略匹配结果，直接把预测位置作为控制点的初始位置。这样既避免了错误匹配造成的影响，又达到了在简单、计算量小条件下较好的预测效果。门限要根据匹配块大小，以及实际图像情况进行调整，图像变化以及噪声剧烈的，门限就要高一点，变化及噪声平稳的就小一点。通常取值为10-20。通过使用本发明方法可以达到图6所示的效果并可以与图3、图4、图5进行比较。