CN102592287B - 基于3d视频的时空域运动分割与估计模型的凸优化方法 - Google Patents

Abstract

本发明公开了一种基于3D视频的时空域运动分割与估计模型的凸优化方法，实施步骤如下：1)根据活动轮廓思想以及背景三维运动参数与二维光流的映射关系建立基于3D视频的时空域运动分割与估计模型；2)将模型转换为对应的水平集描述方程，求出对应的梯度下降方程，求出所述梯度下降方程的等价方程，求出对应于所述等价方程的能量泛函，对能量泛函进行凸松弛，得到凸优化后的时空域运动分割与估计模型；3)引入代价变量对凸优化后的时空域运动分割与估计模型进行进一步的松弛，采用多变量交替迭代算法最小化凸优化后的时空域运动分割与估计模型，迭代收敛后根据选择的门限函数得到最终的分割曲面。本发明具有对目标数量变化的适应性强、分割结果不依赖初始化轮廓、运算效率高的优点。

Description

基于3D视频的时空域运动分割与估计模型的凸优化方法

技术领域

本发明涉及基于计算机视觉的运动分析领域，特别涉及一种利用凸松弛的概念将基于3D视频的运动分割与估计模型转化为全局凸优化极值问题的方法。

背景技术

视频序列的运动分析是计算机视觉领域中的一个基本研究课题，它的应用领域包括治安监控、机器视觉、自动导航、国防兵器、数字媒体、视频编码、3D电视、虚拟现实和智能交通等领域，这些都是国内外重点发展和研究领域。

尽管从二维序列图像或视频中进行运动分割、提取场景获得三维结构和物体的运动信息，无论从理论上还是在实际中都很重要，但由于图像形成的复杂性、三维运动场与图像特征的二维运动场之间的非线性关系和缺乏深度信息不能够对二维运动场精确建模等原因，所得的解通常是在某种先验知识作为约束条件下的近似解，且所得解对数据的噪声非常敏感。随着立体视觉系统和距离传感器的发展和进入实用，三维数据能够比较容易的直接获得。由于三维数据的信息量远高于二维图像，能够准确地表达三维运动和二维运动，三维运动的估计大为简化，将三维数据用于运动分割、提取场景的三维结构和物体的运动信息估计，可以期望得到更为准确和鲁棒的解。

Osher和Sethian首先提出了依赖于时间的运动曲面的水平集描述。水平集方法避免了对拓扑结构变化的处理，计算稳定，已在图像/视频处理和计算机视觉等领域得到了广泛的应用。国内外，不少学者将水平集方法应用在运动分割和光流估算上。有文献将由摄像机运动引起的背景运动用匀速模型表示，同时实现由摄像机运动引起的二维运动场估算和分割曲面的演化，曲面的演变用水平集方法实现。我们在前面研究了基于偏微分方程组水平集方法的多运动目标分割、跟踪、背景运动补偿和运动速度估计，方法中考虑了背景各点运动速度不相同的情况，它能同时进行背景光流的估算、背景运动补偿和运动目标的分割，但不能对摄像机的三维运动进行估算。然而，水平集方法也有一定的局限性，就是它对初始化条件敏感，产生这种现象的原因正是因为能量函数的非凸性使得水平集表示存在局部最小值。目前，很多的研究者将目光放在了如何建立全局凸模型的方法上，克服了水平集方法局部最小值和速度慢的问题，成为又一个研究热点。

发明内容

本发明要解决的技术问题是提供一种对目标数量变化的适应性强、分割结果不依赖初始化轮廓、运算效率高的基于3D视频的时空域运动分割与估计模型的凸优化方法。

为解决上述技术问题，本发明采用的技术方案为：

一种基于3D视频的时空域运动分割与估计模型的凸优化方法，其实施步骤如下：

1)根据活动轮廓思想以及背景三维运动参数与二维光流的映射关系建立基于3D视频的时空域运动分割与估计模型；

2)将所述时空域运动分割与估计模型转换为对应的水平集描述方程，求出所述水平集描述方程对应的梯度下降方程，求出所述梯度下降方程的等价方程，求解对应于所述等价方程的能量泛函，对能量泛函进行凸松弛，得到凸优化后的时空域运动分割与估计模型；

3)引入代价变量对凸优化后的时空域运动分割与估计模型进行进一步的松弛，采用多变量交替迭代算法最小化所述凸优化后的时空域运动分割与估计模型，迭代收敛后根据选择的门限函数得到最终的分割曲面。

作为本发明上述技术方案的进一步改进：

所述步骤1)中建立的时空域运动分割与估计模型的函数表达式为：

$E=\underset{C}{\∫}\mathrm{ds}+\mathrm{\λ}\underset{C^{\mathrm{in}}}{\∫}{e}^{{-(I_t+s\·T+q\·W)}^2}\mathrm{dx}+\mathrm{\λ}\underset{C^{\mathrm{out}}}{\∫}\mathrm{\α}{(I_t+s\·T+q\·W)}^2\mathrm{dx},$

其中E为能量泛函，C为分割曲面，T为背景的平移变量，W为背景的旋转变量，T和W的表达分别为T＝(t₁，t₂，t₃)、W＝(w₁，w₂，w₃)，I_t为图像的时间差分，Cⁱⁿ、C^out分别为时空域中曲面的内部(目标)和外部(背景)区域，ds为单位曲面，dx为单位像素点，α、λ为比例因子，s和q分别为对应于每个像素点的常数向量；s和q的函数表达式分别为：

$s=\left(\begin{array}{c}\frac{{\mathrm{fI}}_x}{Z}\\ \frac{{\mathrm{fI}}_y}{Z}\\ \frac{-{\mathrm{xI}}_x-y{I}_y}{Z}\end{array}\right),$ $q=\left(\begin{array}{c}-{\mathrm{fI}}_y-\frac{y}{f}({\mathrm{xI}}_x+{\mathrm{yI}}_y)\\ {\mathrm{fI}}_x+\frac{x}{f}({\mathrm{xI}}_x+{\mathrm{yI}}_y)\\ -{\mathrm{yI}}_x+{\mathrm{xI}}_y\\ \end{array}\right)$

其中，I_x为图像的横向差分，I_y为图像的纵向差分，x为像素点在图像坐标中的x轴坐标值，y为像素点在图像坐标中的y轴坐标值，Z为像素点在空间坐标中的z轴坐标值，f为摄像机焦距。

所述步骤2)的详细步骤为：

A)将所述时空域运动分割与估计模型转换为对应的水平集描述方程，所述水平集描述方程的函数表达式为：

$\underset{\mathrm{\φ},(T,W)}{\min }\{E_{\mathrm{LSM}}=\underset{\mathrm{\Ω}}{\∫}|\▿\mathrm{\φ}|\mathrm{\δ}\left(\mathrm{\φ}\right)+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{in}}(T,W,x)H\left(\mathrm{\φ}\right)+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{out}}(T,W,x)(1-H\left(\mathrm{\φ}\right))\}$

其中，Ω为整个时空域积分区域，

为内部(目标)区域的运动描述符，

为外部(背景)区域的运动描述符，φ为符号距离函数，

为符号距离函数φ的梯度的模，δ(φ)为关于符号距离函数φ的Dirac函数，H(φ)为关于符号距离函数φ的Heaviside函数，λ为比例因子；

和

的函数表达式为：

$g_r^{\mathrm{in}}(T,W,x)=e^{-{(I_t+s\·T+q\·W)}^2},$ $g_r^{\mathrm{out}}(T,W,x)=\mathrm{\α}{(I_t+s\·T+q\·W)}^2,$

其中，I_t为图像的时间差分，T为背景的平移变量，W为背景的旋转变量，α为比例因子，s和q分别为对应于每个像素点的常数向量；

B)求出所述水平集描述方程对应的梯度下降方程，所述梯度下降方程的函数表达式为：

$\frac{\∂\mathrm{\φ}}{\∂t}=(\mathrm{div}\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}+\mathrm{\λ}(g_r^{\mathrm{in}}(T,W,x)-g_r^{\mathrm{out}}(T,W,x)))\mathrm{\δ}\left(\mathrm{\φ}\right)$

其中T为背景的平移变量，W为背景的旋转变量，x为图像空间中某一个像素点x，t为时间步长，为符号距离函数φ对时间步长的偏微分，

为符号距离函数φ的散度，λ为比例因子，

为内部区域的运动描述符，

为外部区域的运动描述符，δ(φ)为关于符号距离函数φ的Dirac函数；

C)由于δ≥0，求出所述梯度下降方程的等价方程：

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{div}\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}+\mathrm{\λ}(g_r^{\mathrm{in}}(T,W,x)-g_r^{\mathrm{out}}(T,W,x))$

求解对应于所述梯度下降方程的等价方程的能量泛函：

$\underset{\mathrm{\φ}\∈\left\{\mathrm{0,1}\right\}}{\min }\{E=\underset{\mathrm{\Ω}}{\∫}|\▿\mathrm{\φ}|\mathrm{dx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{in}}\mathrm{\φdx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{out}}(1-\mathrm{\φ})\mathrm{dx}\}$

将φ∈{0，1}松弛到连续区间[0，1]上，并用标记符号u代替φ，得到凸优化后的时空域运动分割与估计模型，所述凸优化后的时空域运动分割与估计模型的函数表达式为：

$\underset{u\∈\left[\mathrm{0,1}\right]}{\min }\{E_G=\underset{\mathrm{\Ω}}{\∫}|\▿u|\mathrm{dx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}r(T,W,x)\mathrm{udx}\}$

其中，E_G为凸优化后的能量泛函，u为分割变量，

为分割变量的梯度模值，r(T，W，x)为区域运动描述符，λ为比例因子，x为图像空间中某一个像素点x；r(T，W，x)的函数表达式为 $r(T,W,x)=g_r^{\mathrm{in}}-g_r^{\mathrm{out}}=e^{-{(I_t+s\·T+q\·W)}^2}-\mathrm{\α}{(I_t+s\·T+q\·W)}^2.$

所述步骤3)的详细步骤为：

a)引入一个辅助变量v∈[0，1]来松弛对分割变量的限制，将所述优化后的时空域运动分割与估计模型进一步松弛，进一步松弛后的时空域运动分割与估计模型的函数表达式为：

$E_G(u,v,(T,W))=\underset{\mathrm{\Ω}}{\∫}\left(\right|\▿u|+\mathrm{\λr}(T,W,x)v+\frac{1}{2\mathrm{\ϵ}}{\left|\right|u-v\left|\right|}_{L2}^2)\mathrm{dx};$

b)将背景的平移变量T和旋转变量W初始化为0，将分割变量u和代价变量v初始化为任意值；

c)固定分割变量和代价变量，更新运动变量和区域运动描述符；

d)引入一个与分割变量对应的对偶变量，固定运动变量和代价变量，利用原始-对偶变量的对应关系更新分割变量；

e)固定运动变量和分割变量，更新代价变量；然后重复迭代所述步骤c)～步骤e)直至所述进一步松弛后的时空域运动分割与估计模型收敛；

f)选择门限函数σ∈(0，1)，根据所述门限函数σ∈(0，1)得到最终的分割曲面C＝{x|u(x)＝σ}。

所述步骤c)的详细步骤为：

对背景的运动参数(T，W)建立超定线性方程组θ_iA_i＝b_i，其中θ_i＝(T_i，W_i)为第i次迭代的背景运动变量，A_i是6×N_i的矩阵，b_i是1×N_i的向量，其中N_i为第i次迭代时背景区域的像素点个数， $A_i=(\left[\begin{array}{c}s\left(1\right)\\ q\left(1\right)\end{array}\right],...,\left[\begin{array}{c}s\left(N_i\right)\\ q\left(N_i\right)\end{array}\right]),$ b_i＝(-I_t(1)，...，-I_t(N_i))，然后采用最小二乘法求解所述超定线性方程组得到运动参数(T，W)的估计值，并根据所述估计值更新

和r(T，W，x)；

所述步骤d)的详细步骤为：

引入与原始分割变量u相对应的对偶变量

当运动变量和代价变量固定，根据原始-对偶变量定义，u的泛函极值子问题的函数表达式为：

求解所述泛函极值子问题得到u的迭代公式

其中uⁱ⁺¹为u在第i+1次迭代的值，vⁱ为v在第i次迭代的值，

为梯度算子，

为

在第i+1次迭代的值，为梯度算子，ε为常数变量；

利用半隐式梯度下降法，得到对偶变量

的迭代式：

其中，

v为代价变量，δt为时间步长。

所述步骤e)的详细步骤为：

当运动变量和分割变量固定时，代价变量的泛函极值子问题为：

$v=\arg \underset{v}{\min }\underset{\mathrm{\Ω}}{\∫}(\mathrm{\λr}(T,W,x)v+\frac{1}{2\mathrm{\ϵ}}{\left|\right|u-v\left|\right|}_{L2}^2)\mathrm{dx}$

对所述极值子问题求解，得到代价变量的迭代式为：

vⁱ⁺¹＝min{max{uⁱ⁺¹-ελrⁱ⁺¹(T，W，x)，0}，1}

其中vⁱ⁺¹为v在第i+1次迭代的值，uⁱ⁺¹为u在第i+1次迭代的值，ε为常数变量，λ为比例因子，rⁱ⁺¹(T，W，x)为区域运动描述符在第i+1次迭代的值。

本发明具有下述优点：

1、本发明通过活动轮廓思想以及背景三维运动参数与二维光流的映射关系建立基于3D视频的时空域运动分割与估计模型，能够实现基于时空域的动态背景情况下的多目标跟踪，能有效的利用空间和时域信息更好施加三维运动约束，适合于处理目标数目有变化的情况。

2、本发明通过对能量泛函进行凸松弛得到凸优化后的时空域运动分割与估计模型，然后通过引入代价变量对凸优化后的时空域运动分割与估计模型进行进一步的松弛，实现对时空域运动分割与估计模型的全局凸优化，能够避免初始化轮廓对结果的影响。

3、本发明通过引入代价变量对凸优化后的时空域运动分割与估计模型进行进一步的松弛，并引入对偶变量，采用多变量交替迭代算法最小化凸优化后的时空域运动分割与估计模型，极大地提高了运算效率。

附图说明

图1为本发明实施例的基本流程示意图。

图2为本发明实施例步骤3)的基本流程示意图。

图3为本发明实施例中3D视频在时空域运动分割与估计模型收敛后对某一帧图像估计出的光流场示意图。

图4为本发明实施例中基于图3的图像帧的分割变量演化结果示意图。

图5为本发明实施例中基于图3的图像帧的最终得到的分割曲面示意图。

具体实施方式

如图1所示，本发明实施例基于3D视频的时空域运动分割与估计模型的凸优化方法的实施步骤如下：

1)根据活动轮廓思想以及背景三维运动参数与二维光流的映射关系建立基于3D视频的时空域运动分割与估计模型；

2)将时空域运动分割与估计模型转换为对应的水平集描述方程，求出水平集描述方程对应的梯度下降方程，求出梯度下降方程的等价方程，求解对应于等价方程的能量泛函，对能量泛函进行凸松弛，得到凸优化后的时空域运动分割与估计模型；

3)引入代价变量对凸优化后的时空域运动分割与估计模型进行进一步的松弛，采用多变量交替迭代算法最小化凸优化后的时空域运动分割与估计模型，迭代收敛后根据选择的门限函数得到最终的分割曲面。

本实施例中，步骤1)中建立的时空域运动分割与估计模型的函数表达式为：

$E=\underset{C}{\∫}\mathrm{ds}+\mathrm{\λ}\underset{C^{\mathrm{in}}}{\∫}{e}^{{-(I_t+s\·T+q\·W)}^2}\mathrm{dx}+\mathrm{\λ}\underset{C^{\mathrm{out}}}{\∫}\mathrm{\α}{(I_t+s\·T+q\·W)}^2\mathrm{dx},$

其中E为能量泛函，C为分割曲面，T为背景的平移变量，W为背景的旋转变量，T和W的表达分别为T＝(t₁，t₂，t₃)、W＝(w₁，w₂，w₃)，I_t为图像的时间差分，Cⁱⁿ、C^out分别为时空域中曲面的内部(目标)和外部(背景)区域，ds为单位曲面，dx为单位像素点，α、λ为比例因子，s和q分别为对应于每个像素点的常数向量；s和q的函数表达式分别为：

$s=\left(\begin{array}{c}\frac{{\mathrm{fI}}_x}{Z}\\ \frac{{\mathrm{fI}}_y}{Z}\\ \frac{-{\mathrm{xI}}_x-y{I}_y}{Z}\end{array}\right),$ $q=\left(\begin{array}{c}-{\mathrm{fI}}_y-\frac{y}{f}({\mathrm{xI}}_x+{\mathrm{yI}}_y)\\ {\mathrm{fI}}_x+\frac{x}{f}({\mathrm{xI}}_x+{\mathrm{yI}}_y)\\ -{\mathrm{yI}}_x+{\mathrm{xI}}_y\\ \end{array}\right)$

其中，I_x为图像的横向差分，I_y为图像的纵向差分，x为像素点在图像坐标中的x轴坐标值，y为像素点在图像坐标中的y轴坐标值，Z为像素点在空间坐标中的z轴坐标值，f为摄像机焦距。

一般来说，两相的主动轮廓模型可以表示为：

$\min \{E_{\mathrm{AC}}\left(C\right)=\underset{C}{\∫}{g}_b(C,s)\mathrm{ds}+\mathrm{\λ}\underset{C^{\mathrm{in}}}{\∫}{g}_r^{\mathrm{in}}(C^{\mathrm{in}},x)\mathrm{dx}+\mathrm{\λ}\underset{C^{\mathrm{out}}}{\∫}{g}_r^{\mathrm{out}}(C^{\mathrm{out}},x)\mathrm{dx}\}$ (1)

其中C代表闭合曲线(2D图像)或曲面(3D图像)，Cⁱⁿ和C^out分别代表图像空间中C的内部和外部区域，g_b：Ω→R是边缘函数(比如测地线模型中的边缘检测函数)，

分别为自定义的内部和外部区域函数，ds、dx分别为单位长度/面积和单位像素点，λ为权重参数。

在摄像机运动的情况下建立对目标的运动分割与估计模型，首先是选择适当的参数模型对背景运动建模。假设背景为一个运动刚体，其三维运动(T，W)分别用平移变量T＝(t₁，t₂，t₃)和旋转变量W＝(w₁，w₂，w₃)表示。背景中任意一点运动产生的光流w＝(u，v)满足光流约束方程

其中

代表图像梯度。根据光流场与三维运动参数的映射关系，背景区域的任意像素点满足以下三维光流约束方程：

I_op(T，W，x)＝I_t+s·T+q·W＝0    (2)

其中s和q均为3*1的常向量：

$s=\left(\begin{array}{c}\frac{{\mathrm{fI}}_x}{Z}\\ \frac{{\mathrm{fI}}_y}{Z}\\ \frac{-{\mathrm{xI}}_x-y{I}_y}{Z}\end{array}\right),$ $q=\left(\begin{array}{c}-{\mathrm{fI}}_y-\frac{y}{f}({\mathrm{xI}}_x+{\mathrm{yI}}_y)\\ {\mathrm{fI}}_x+\frac{x}{f}({\mathrm{xI}}_x+{\mathrm{yI}}_y)\\ -{\mathrm{yI}}_x+{\mathrm{xI}}_y\\ \end{array}\right)$

当(T，W)正确估计时，背景区域的I_op应该趋近于0；反之，前景区域的I_op就较大，因此3D光流约束方程是一个很好的用于区别前景和背景的特征。我们可以选用以下观察模型来设计区域运动函数：

$\begin{array}{cc}{g}_r^{\mathrm{in}}(T,W,x)=e^{-I_{\mathrm{op}}^2(T,W,x)}& x\∈C^{\mathrm{in}}\end{array}$ (3)

$\begin{array}{cc}{g}_r^{\mathrm{out}}(T,W,x)=\mathrm{\α}{I}_{\mathrm{op}}^2(T,W,x)& x\∈C^{\mathrm{out}}\end{array}$

其中α是用于控制内部和外部区间竞争的变量。另外，g_b＝1定义为时空曲面上的惩罚项，用来平滑三维曲面。这样，我们就得到了基于3D序列的时空域运动分割与估计模型：

$E=\underset{C}{\∫}\mathrm{ds}+\mathrm{\λ}\underset{C^{\mathrm{in}}}{\∫}{e}^{{-(I_t+s\·T+q\·W)}^2}\mathrm{dx}+\mathrm{\λ}\underset{C^{\mathrm{out}}}{\∫}\mathrm{\α}{(I_t+s\·T+q\·W)}^2\mathrm{dx}$ (4)

本实施例中，步骤2)的详细步骤为：

A)将时空域运动分割与估计模型转换为对应的水平集描述方程，水平集描述方程的函数表达式为：

$\underset{\mathrm{\φ},(T,W)}{\min }\{E_{\mathrm{LSM}}=\underset{\mathrm{\Ω}}{\∫}|\▿\mathrm{\φ}|\mathrm{\δ}\left(\mathrm{\φ}\right)+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{in}}(T,W,x)H\left(\mathrm{\φ}\right)+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{out}}(T,W,x)(1-H\left(\mathrm{\φ}\right))\}$

其中，Ω为整个时空域积分区域，

为内部(目标)区域的运动描述符，

为外部(背景)区域的运动描述符，φ为符号距离函数，

为符号距离函数φ的梯度的模，δ(φ)为关于符号距离函数φ的Dirac函数，H(φ)为关于符号距离函数φ的Heaviside函数，λ为比例因子；

和

的函数表达式为：

$g_r^{\mathrm{in}}(T,W,x)=e^{-{(I_t+s\·T+q\·W)}^2},$ $g_r^{\mathrm{out}}(T,W,x)=\mathrm{\α}{(I_t+s\·T+q\·W)}^2,$

其中，I_t为图像的时间差分，T为背景的平移变量，W为背景的旋转变量，α为比例因子，s和q分别为对应于每个像素点的常数向量；

B)求出水平集描述方程对应的梯度下降方程，梯度下降方程的函数表达式为：

$\frac{\∂\mathrm{\φ}}{\∂t}=(\mathrm{div}\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}+\mathrm{\λ}(g_r^{\mathrm{in}}(T,W,x)-g_r^{\mathrm{out}}(T,W,x)))\mathrm{\δ}\left(\mathrm{\φ}\right)$

其中T为背景的平移变量，W为背景的旋转变量，x为图像空间中某一个像素点x，t为时间步长，为符号距离函数φ对时间步长的偏微分，

为符号距离函数φ的散度，λ为比例因子，

为内部区域的运动描述符，

为外部区域的运动描述符，δ(φ)为关于符号距离函数φ的Dirac函数；

C)由于δ≥0，求出梯度下降方程的等价方程：

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{div}\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}+\mathrm{\λ}(g_r^{\mathrm{in}}(T,W,x)-g_r^{\mathrm{out}}(T,W,x))$

求解对应于梯度下降方程的等价方程的能量泛函：

$\underset{\mathrm{\φ}\∈\left\{\mathrm{0,1}\right\}}{\min }\{E=\underset{\mathrm{\Ω}}{\∫}|\▿\mathrm{\φ}|\mathrm{dx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{in}}\mathrm{\φdx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{out}}(1-\mathrm{\φ})\mathrm{dx}\}$

将φ∈{0，1}松弛到连续区间[0，1]上，并用标记符号u代替φ，得到凸优化后的时空域运动分割与估计模型，凸优化后的时空域运动分割与估计模型的函数表达式为：

$\underset{u\∈\left[\mathrm{0,1}\right]}{\min }\{E_G=\underset{\mathrm{\Ω}}{\∫}|\▿u|\mathrm{dx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}r(T,W,x)\mathrm{udx}\}$

其中，E_G为凸优化后的能量泛函，u为分割变量，

为分割变量的梯度模值，r(T，W，x)为区域运动描述符，λ为比例因子，x为图像空间中某一个像素点x；r(T，W，x)的函数表达式为 $r(T,W,x)=g_r^{\mathrm{in}}-g_r^{\mathrm{out}}=e^{-{(I_t+s\·T+q\·W)}^2}-\mathrm{\α}{(I_t+s\·T+q\·W)}^2.$

能量函数的最小化通常使用水平集方法进行求解，它的主要思想是用n+1维函数的零水平集隐性的表示n维曲线C：Ω→Rⁿ，这种方法的优势是可以处理变化的拓扑结构比如多条曲线的分裂和合并，且计算稳定。然而，水平集方法也有以下不足之处：1)每次迭代后都需要重新初始化符号距离函数，计算效率低；2)因为存在局部最小值问题，所有对初始化曲线敏感。为了克服上述问题，我们将对得到的3D时空域运动分割与估计模型进行凸优化。

基于水平集方法的3D时空域运动分割与估计模型(4)的目标函数可以表示为：

$\underset{\mathrm{\φ},(T,W)}{\min }\{E_{\mathrm{LSM}}=\underset{\mathrm{\Ω}}{\∫}|\▿\mathrm{\φ}|\mathrm{\δ}\left(\mathrm{\φ}\right)+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{in}}(T,W,x)H\left(\mathrm{\φ}\right)+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{out}}(T,W,x)(1-H\left(\mathrm{\φ}\right))\}$ (5)

其中Ω为整个时空域积分区域，φ为以符号距离函数表示的水平集函数，H(φ)和δ(φ)分别为Heaviside函数和Dirac函数，和

为前面定义的区域函数。形如式(5)的多变量优化问题通常采用子极值问题的交替优化来实现，当更新当前的后，水平集函数对应的梯度下降方程为：

$\frac{\∂\mathrm{\φ}}{\∂t}=(\mathrm{div}\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}+\mathrm{\λ}(g_r^{\mathrm{in}}(T,W,x)-g_r^{\mathrm{out}}(T,W,x)))\mathrm{\δ}\left(\mathrm{\φ}\right)$ (6)

由于δ≥0且不会影响水平集演化方向，式(6)可以改写成以下等价形式：

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{div}\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}+\mathrm{\λ}(g_r^{\mathrm{in}}(T,W,x)-g_r^{\mathrm{out}}(T,W,x))$ (7)

式(7)正是如下能量泛函对应的梯度下降方程

$\underset{\mathrm{\φ}\∈\left\{\mathrm{0,1}\right\}}{\min }\{E=\underset{\mathrm{\Ω}}{\∫}|\▿\mathrm{\φ}|\mathrm{dx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{in}}\mathrm{\φdx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{g}_r^{\mathrm{out}}(1-\mathrm{\φ})\mathrm{dx}\}$ (8)

在水平集描述中φ的定义区间为{0，1}，由于能量函数(8)定义在非凸集上，故其仍为非全局最优问题。将φ∈{0，1}松弛到连续区间[0，1]上，为了避免混淆，用符号u代替φ，上述问题就转化为如下全局最优模型：

$\underset{u\∈\left[\mathrm{0,1}\right]}{\min }\{E_G=\underset{\mathrm{\Ω}}{\∫}|\▿u|\mathrm{dx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}r(T,W,x)\mathrm{udx}\}$ (9)

其中 $r(T,W,x)=g_r^{\mathrm{in}}-g_r^{\mathrm{out}}=e^{-{(I_t+s\·T+q\·W)}^2}-\mathrm{\α}{(I_t+s\·T+q\·W)}^2.$

如图2所示，本实施例步骤3)的详细步骤为：

a)引入一个辅助变量v∈[0，1]来松弛对分割变量的限制，将优化后的时空域运动分割与估计模型进一步松弛，进一步松弛后的时空域运动分割与估计模型的函数表达式为：

$E_G(u,v,(T,W))=\underset{\mathrm{\Ω}}{\∫}\left(\right|\▿u|+\mathrm{\λr}(T,W,x)v+\frac{1}{2\mathrm{\ϵ}}{\left|\right|u-v\left|\right|}_{L2}^2)\mathrm{dx}.$

b)将背景的平移变量T和旋转变量W初始化为0，将分割变量u和代价变量v初始化为任意值。

c)固定分割变量和代价变量，更新运动变量和区域运动描述符。

d)引入一个与分割变量对应的对偶变量，固定运动变量和代价变量，利用原始-对偶变量的对应关系更新分割变量。

e)固定运动变量和分割变量，更新代价变量；然后重复迭代步骤c)～步骤e)直至进一步松弛后的时空域运动分割与估计模型收敛；本实施例中3D视频中某一帧图像在时空域运动分割与估计模型收敛后估计出的光流场如图3所示，该帧图像的分割变量演化结果图如图4所示。

f)选择门限函数σ∈(0，1)，根据门限函数σ∈(0，1)以及C_in＝{x|u(x)＞σ}，得到最终的分割曲面(如图5所示)C＝{x|u(x)＝σ}。

本实施例中，步骤c)的详细步骤为：

对背景的运动参数(T，W)建立超定线性方程组θ_iA_i＝b_i，其中θ_i＝(T_i，W_i)为第i次迭代的背景运动变量，A_i是6×N_i的矩阵，b_i是1×N_i的向量，其中N_i为第i次迭代时背景区域的像素点个数， $A_i=(\left[\begin{array}{c}s\left(1\right)\\ q\left(1\right)\end{array}\right],...,\left[\begin{array}{c}s\left(N_i\right)\\ q\left(N_i\right)\end{array}\right]),$ b_i＝(-I_t(1)，...，-I_t(N_i))，然后采用最小二乘法求解超定线性方程组得到运动参数(T，W)的估计值，并根据估计值更新

和r(T，W，x)；

本实施例中，步骤d)的详细步骤为：

引入与原始分割变量u相对应的对偶变量

当运动变量和代价变量固定，根据原始-对偶变量定义，u的泛函极值子问题的函数表达式为

求解泛函极值子问题得到u的迭代公式

其中uⁱ⁺¹为u在第i+1次迭代的值，vⁱ为v在第i次迭代的值，

为梯度算子，

为

在第i+1次迭代的值，

为梯度算子，ε为常数变量；

利用半隐式梯度下降法，得到对偶变量

的迭代式：

其中，

v为代价变量，δt为时间步长。

为了求解分割变量，我们引入一个变量

则全变分可以改写为以下形式：

(10)

这里u和

就分别被称为原始和对偶变量。当(T，W)和v固定时，u的泛函极值子问题可以写成：

(11)

根据对偶矢量的定义，对方程(11)求解可得：

(12)

将式(12)代入方程(11)，得到：

(13)

式(13)等同于以下最小化问题：

(14)

根据欧拉-拉格朗日方程，可以得到x点处的最优解的必要条件：

(15)

这里α(x)是满足限制条件

的拉格朗日乘子：

(16)

将(16)代入(15)，得到：

(17)

然后利用半隐式梯度下降法，最后得到对偶变量

的迭代式：

(19)

其中时间步长取δt＝1/8。

代价变量的作用是对足够小的ε，使得u≈v。本实施例中，步骤e)的详细步骤为：

当运动变量和分割变量固定时，代价变量的泛函极值子问题为：

$v=\arg \underset{v}{\min }\underset{\mathrm{\Ω}}{\∫}(\mathrm{\λr}(T,W,x)v+\frac{1}{2\mathrm{\ϵ}}{\left|\right|u-v\left|\right|}_{L2}^2)\mathrm{dx}$

对极值子问题求解，得到代价变量的迭代式为：

vⁱ⁺¹＝min{max{uⁱ⁺¹-ελrⁱ⁺¹(T，W，x)，0}，1}

其中vⁱ⁺¹为v在第i+1次迭代的值，uⁱ⁺¹为u在第i+1次迭代的值，ε为常数变量，λ为比例因子，rⁱ⁺¹(T，W，x)为区域运动描述符在第i+1次迭代的值。

以上所述仅为本发明的优选实施方式，本发明的保护范围并不仅限于上述实施方式，凡是属于本发明原理的技术方案均属于本发明的保护范围。对于本领域的技术人员而言，在不脱离本发明的原理的前提下进行的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种基于3D视频的时空域运动分割与估计模型的凸优化方法，其特征在于其实施步骤如下： 

1）根据活动轮廓思想以及背景三维运动参数与二维光流的映射关系建立基于3D视频的时空域运动分割与估计模型； 

所述的时空域运动分割与估计模型的函数表达式为： 

其中E为能量泛函，C为分割曲面，T为背景的平移变量，W为背景的旋转变量，T和W的表达分别为T＝(t₁,t₂,t₃)、W＝(w₁,w₂,w₃)，I_t为图像的时间差分，Cⁱⁿ、C^out分别为时空域中曲面的内部和外部区域，ds为单位曲面，dx为单位像素点，α、λ为比例因子，s和q分别为对应于每个像素点的常数向量；s和q的函数表达式分别为： 

其中，I_x为图像的横向差分，I_y为图像的纵向差分，x为像素点在图像坐标中的x轴坐标值，y为像素点在图像坐标中的y轴坐标值，Z为像素点在空间坐标中的z轴坐标值，f为摄像机焦距； 

2）将所述时空域运动分割与估计模型转换为对应的水平集描述方程，求出所述水平集描述方程对应的梯度下降方程，求出所述梯度下降方程的等价方程，求解对应于所述等价方程的能量泛函，对能量泛函进行凸松弛，得到凸优化后的时空域运动分割与估计模型； 

3）引入代价变量对凸优化后的时空域运动分割与估计模型进行进一步的松弛，采用多变量交替迭代算法最小化所述凸优化后的时空域运动分割与估计模型，迭代收敛后根据选择的门限函数得到最终的分割曲面。 

2.根据权利要求1所述的基于3D视频的时空域运动分割与估计模型的凸优化方法，其特征在于，所述步骤2）的详细步骤为： 

A）将所述时空域运动分割与估计模型转换为对应的水平集描述方程，所述水平集描述方程的函数表达式为： 

其中，Ω为整个时空域积分区域，

为内部区域的运动描述符，

为外部区域的运动描述符，φ为符号距离函数，|▽φ|为符号距离函数φ的梯度的模，δ(φ)为关于符号距离函数φ的Dirac函数，H(φ)为关于符号距离函数φ的Heaviside函数，λ为比例因子；

和

的函数表达式为： 

其中，I_t为图像的时间差分，T为背景的平移变量，W为背景的旋转变量，α为比例因子，s和q分别为对应于每个像素点的常数向量； 

B）求出所述水平集描述方程对应的梯度下降方程，所述梯度下降方程的函数表达式为： 

其中T为背景的平移变量，W为背景的旋转变量，x为图像空间中某一个像素点x，t为时间步长，

为符号距离函数φ对时间步长的偏微分，

为符号距离函数φ的散度，λ为比例因子，

为内部区域的运动描述符，

为外部区域的运动描述符，δ(φ)为关于符号距离函数φ的Dirac函数； 

C)由于δ≥0，求出所述梯度下降方程的等价方程： 

求解对应于所述梯度下降方程的等价方程的能量泛函： 

将φ∈{0,1}松弛到连续区间[0,1]上，并用标记符号u代替φ，得到凸优化后的时空域运动分割与估计模型，所述凸优化后的时空域运动分割与估计模型的函数表达式为： 

其中，E_G为凸优化后的能量泛函，u为分割变量，|▽u|为分割变量的梯度模值，r(T,W,x)为区域运动描述符，λ为比例因子，x为图像空间中某一个像素点x；r(T,W,x)的函数表达式为

3.根据权利要求2所述的基于3D视频的时空域运动分割与估计模型的凸优化方法，其特征在于，所述步骤3）的详细步骤为： 

a）引入一个代价变量v∈[0,1]来松弛对分割变量的限制，将所述优化后的时空域运动分割与估计模型进一步松弛，进一步松弛后的时空域运动分割与估计模型的函数表达式为： 

b）将背景的平移变量T和旋转变量W初始化为0，将分割变量u和代价变量v初始化为任意值； 

c）固定分割变量和代价变量，更新运动变量和区域运动描述符，其中第i次迭代时背景的运动变量定义为θ_i＝(T_i,W_i)，T_i为第i次迭代时背景的平移变量，W_i为第i次迭代时背景的旋转变量； 

d）引入一个与分割变量对应的对偶变量，固定运动变量和代价变量，利用原始分割变量-对偶变量的对应关系更新分割变量； 

e）固定运动变量和分割变量，更新代价变量；然后重复迭代所述步骤c）～步骤e）直至所述进一步松弛后的时空域运动分割与估计模型收敛； 

f）选择门限函数σ∈(0,1)，根据所述门限函数σ∈(0,1)得到最终的分割曲面C＝{x|u(x)＝σ}。 

4.根据权利要求3所述的基于3D视频的时空域运动分割与估计模型的凸优化方法，其特征在于，所述步骤c）的详细步骤为： 

对背景的运动参数(T,W)建立超定线性方程组θ_iA_i＝b_i，A_i是6×N_i的矩阵，b_i是1×N_i的向量，其中N_i为第i次迭代时背景区域的像素点个数，b_i＝(-I_t(1),...,-I_t(N_i))，然后采用最小二乘法求解所述超定线性方程组得到运动参数(T,W)的估计值，并根据所述估计值更新

和

5.根据权利要求3所述的基于3D视频的时空域运动分割与估计模型的凸优化方法，其特征在于，所述步骤d）的详细步骤为： 

引入与原始分割变量u相对应的对偶变量

当运动变量和代价变量固定，根据原始-对偶变量定义，u的泛函极值子问题的函数表达式为： 

求解所述泛函极值子问题得到u的迭代公式

其中uⁱ⁺¹为u在第i+1次迭代的值，vⁱ为v在第i次迭代的值，▽为梯度算子，

为

在第i次迭代的值，▽为梯度算子，ε为常数变量； 

利用半隐式梯度下降法，得到对偶变量

的迭代式： 

其中，

v为代价变量，δt为时间步长。 

6.根据权利要求3所述的基于3D视频的时空域运动分割与估计模型的凸优化方法，其特征在于，所述步骤e）的详细步骤为： 

当运动变量和分割变量固定时，代价变量的泛函极值子问题为： 

对所述极值子问题求解，得到代价变量的迭代式为： 

vⁱ⁺¹＝min{max{uⁱ⁺¹-ελrⁱ⁺¹(T,W,x),0},1} 

其中vⁱ⁺¹为v在第i+1次迭代的值，uⁱ⁺¹为u在第i+1次迭代的值，ε为常数变量，λ为比例因子，rⁱ⁺¹(T,W,x)为区域运动描述符在第i+1次迭代的值。 

Images