多视环境中基于水平集函数演化的非刚性目标场景流估计

摘要

本发明公开了多视环境中基于水平集函数演化的非刚性目标场景流估计，具体按照以下步骤实施：步骤1：多视环境设置及预处理；步骤2：获取静态目标三维外形；步骤3：构建时间序列上的关于水平集函数的能量泛函；步骤4：获取动态非刚性目标三维外形及其场景流。本发明用水平集函数的演化模拟非刚性目标的运动，在时空颜色相关性约束、投影区域一致性约束下有效地驱动非刚性目标初始外形曲面向目标曲面演化，而水平集函数自身特点的约束，使的水平集函数演化后仍为水平集函数，避免了水平集函数重新初始化的问题，解决了现在估计非刚性目标的场景流时，由于环境因素，计算方法不当，往往会使得场景流不够连续、不平滑或存在空洞的问题。

说明书

技术领域

本发明属于计算机视觉技术领域，具体涉及一种多视环境中基于水平集函数演化的非刚性目标场景流估计。

背景技术

多视环境中非刚性目标的三维外形场景流估计指在多视环境下获取目标在不同视点的图像，依据计算机视觉原理恢复目标的三维外形，然后根据多视环境获取的目标的运动视频，使用计算机视觉、数字视频处理方法估计其三维外形的运动场，即场景流。非刚性目标三维外形场景流估计方法可获取动态非刚性目标不同时刻三维外形的关联性，因此在虚拟现实、人机交互等方面有着重要的应用。

文献Three-Dimensional Scene Flow(IEEE Transactions on Pattern Analysis and Machine Intelligence,27(3):475–480,2005)和Time Varying Volumetric Scene Reconstruction Using Scene Flow(The 18th British Machine Vision Conference,2007)首先使用多视摄像机获取目标的运动视频，然后分别求取每个摄像机获取的视频的光流场，最后依据计算机视觉理论，将多摄像机标定参数与多个摄像机中的光流场结合，恢复目标的场景流。这类方法获取的场景流其实是由不同视点场景流拼合而成的，因此摄像机标定误差、噪声等因素往往会使得最终获取的场景流存在不连续、不平滑现象。

发明内容

本发明的目的是提供一种多视环境中基于水平集函数演化的非刚性目标场景流估计，解决了现有获取的场景流存在不连续、不平滑的问题。

本发明所采用的技术方案是，多视环境中基于水平集函数演化的非刚性目标场景流估计，具体按照以下步骤实施：

步骤1：多视环境设置及预处理；

步骤2：获取静态目标三维外形；

步骤3：构建时间序列上的关于水平集函数的能量泛函；

步骤4：获取动态非刚性目标三维外形及其场景流。

本发明的特点还在于：

步骤1具体为：

将nCamNum个摄像机环绕观测场景放置，并使用张氏标定法获取摄像机的内、外参数；然后使用nCamNum个摄像机采集不包含目标的背景图像{Back₁,Back₂,...,Back_j,...,Back_nCamNum}；之后将动态非刚性目标O置于观测场景中，多个摄像机同步采集目标的运动视频，在时刻t、时刻t+1，同步采集到的视频帧分别为{I_1,t,I_2,t,...,I_j,t,...,I_nCamNum,t}、{I_1,t+1,I_2,t+1,...,I_j,t+1,...,I_nCamNum,t+1}，最后使用背景减法获取目标O的侧影轮廓，即目标图像，在时刻t和时刻t+1，分别使用多视图像帧减去背景图像，从而获取目标图像：{Silhouette_1,t,...,Silhouette_j,t,...,Silhouette_nCamNum,t}、{Silhouette_1,t+1,...,Silhouette_j,t+1,...,Silhouette_nCamNum,t+1}。

步骤2具体为：

将多视摄像机的公共观测区域离散化为一个三维体素空间Ω；

首先，使用{Silhouette_1,t,...,Silhouette_j,t,...,Silhouette_nCamNum,t}，基于体素着色法重建目标三维外形的体素模型Shape_Data_t，其中，

然后，重建目标三维外形，及初始化场景流估计所需的水平集函数：

Step1：初始化水平集函数φ(x)：在Ω内，靠近Ω边缘的区域构造一个足够大的封闭初始曲面Γ(0)，确保其能够包括Shape_Data的所有体素，构造初始水平集函数φ₀(x)：

其中c₀为一预先给定的正整数；

Step2：构造关于水平集函数φ(x)的能量函数E(φ)：

$E\left(\phi \right)=\alpha ·\underset{\Omega }{\int }g\left(x\right)\delta \left(\phi \right)\left|▿\phi \right|dx+\beta ·P\left(\phi \right)+\gamma ·A\left(\phi \right)---\left(2\right)$

其中，α、β和γ为比例系数，g(x)为在Ω中计算距离场为：g(x)＝dist(x,Shape_Data)，其中，x∈Ω，dist(x,Shape_Data)表示体素x到三维外形Shape_Data的最短距离，δ为Dirac函数；

P(φ)为惩罚项为：

$P\left(\phi \right)=\underset{\Omega }{\int }\rho \left(\right|▿\phi \left|\right)dx---\left(3\right)$

为：

$\rho \left(\right|▿\phi \left|\right)=\left(\begin{array}{cc}\frac{1}{{\left(2\pi \right)}^2}·(1-\cos (2\pi ·\left|▿\phi \right|\left)\right)& \left|▿\phi \right|\le 1\\ \frac{1}{2}·{\left(\right|▿\phi |-1)}^2& \left|▿\phi \right|>1\end{array}\right)---\left(4\right)$

A(φ)为加速项为：

$A\left(\phi \right)=\underset{\Omega }{\int }H(g-\eta )g\left(x\right)H(-\phi )dx---\left(5\right)$

H(x)是Heaviside函数，η是偏移程度，当g＜η时，加速项失效，从而可确保曲面不会穿过空洞；

对于能量函数E(φ)，其相应的梯度下降流，即水平集函数的演化速度为：

Step3：令更新水平集函数：

φⁿ⁺¹(x)＝φⁿ(x)+Δτ·D(φⁿ(x))>

n表示迭代次数，Δτ为时间步长；

Step4：判断是否小于阈值来判断是否收敛，若不收敛，则按照式(7)计算φⁿ⁺¹(x)的值，并转至Step4，若收敛，则结束；

令f^t(x)＝φⁿ⁺¹(x)，作为时刻t目标三维外形的水平集函数表示。

在曲面演化开始前，η赋值为以体素为单位的模型最大空洞的直径长度。

为了加快Step3的更新速度，使用窄带法：提取零水平集φ(x,τ)＝0，将提取出的体素放入集合LevelSet中，对于每个体素i，i∈LevelSet，构造窄带其中为体素i的邻域，r为邻域的范围，通常指定其为以i为中心，2r+1为边长的立方体，即更新只在零水平集的邻域内进行。

步骤3具体为：

将多视环境获取的非刚性目标运动视频提供的信息作为约束条件，构建关于水平集函数的能量泛函E(f^t+1,V_t)：

①时空颜色相关性约束项

时间颜色相关性φ_temp(vox，V_t^vox)是指在时刻t目标三维外形上的某一体素以速度V_t^vox运动，在时刻t+1对应的目标三维外形上的体素为vox，其运动前后的颜色保持一致；

时间颜色相关性φ_temp(vox,V_t^vox)为：

$\begin{array}{c}\phi \_temp(vox,V_t^{vox})\\ =\frac{1}{nVisNum}\underset{j=1}{\overset{nVisNum}{\Sigma }}\left(\frac{\partial photo\left({\mathrm{pro}}_j\right(vox),t+1)}{\partial vox}\right)·V_t^{vox}+\frac{\partial photo\left({\mathrm{pro}}_j\right(vox),t+1)}{\partial (t+1)}{)}^2\end{array}---\left(8\right)$

nVisNum表示对体素vox可见的视点个数，pro_j(vox)表示体素vox在视点j的二维投影，photo(X,t+1)为时刻t+1像素X对应的像素颜色值，φ_temp(vox,V_t^vox)反映了目标三维外形上的体素运动前后在多视图像中投影的颜色相关性：φ_temp(vox,V_t^vox)越小，颜色相关性越大；

空间颜色相关性φ_space(vox,t+1)是指时刻t+1目标三维外形上体素vox与其周围的体素具有颜色一致性；

空间颜色相关性φ_space(vox,t+1)为：

$\phi \_space(vox,t+1)=\frac{1}{nVisNum}\underset{j=1}{\overset{nVisNum}{\Sigma }}{\left(photo\right({\mathrm{pro}}_j\left(vox\right),t+1)-\overline{color})}^2---\left(9\right)$

其中：

$\overline{color}=\frac{1}{nVisNum}\underset{j=1}{\overset{nVisNum}{\Sigma }}photo\left({\mathrm{pro}}_j\right(vox),t+1)$

将时间颜色相关性和空间颜色相关性结合，得出时空颜色相关性约束项：

$\begin{array}{c}\phi \_color(vox,V_t^{vox})\\ =isColored(vox,V_t^{vox})·\phi \_temp(vox,V_t^{vox})+(1-isColored(vox,V_t^{vox}\left)\right)·\phi \_space(vox,t+1)\end{array}---\left(10\right)$

其中，isColored(vox,V_t^vox)表示时刻t体素vox-V_t^vox是否被判断为目标三维外形上的体素且已被着色，其取值方式为：

此时，Ω中所有体素的时空颜色相关性约束对应的能量函数为：

$E_{color}(f^{t+1},V_t)=\underset{vox\in \Omega }{\int \int \int }\phi \_color(vox,V_t^{vox})·\left|▿H\right|dxdydz---\left(11\right)$

②投影区域一致性约束项

投影区域一致性约束项对应的能量函数为：

$E_{region}\left(f^{t+1}\right)=\underset{vox\in \Omega }{\int \int \int }(in\_Silhouette(vox)·H(f^{t+1})+(1-in\_Silhouette\left(vox\right))·H(-f^{t+1}\left)\right)dxdydz---\left(12\right)$

其中in_Silhouette(vox)表示体素vox在多视图像中的投影是否都落入目标的侧影轮廓中：

③水平集函数规整项

水平集函数规整项对应的能量函数P(f)为：

$P\left(f^{t+1}\right)=\underset{\Omega }{\int }\rho \left(\right|▿f^{t+1}\left|\right)dx---\left(13\right)$

在上述约束项和规整项的基础上，能量泛函E(f^t+1,V_t)为：

E(f^t+1,V_t)＝κ·E_color(f^t+1,V_t)+λ·E_region(f^t+1)+μ·P(f^t+1)>

其中，μ＞0为常量参数，用于控制f偏离符号距离函数后的惩罚效应，κ、λ表示权重系数。

步骤4具体为：

最小化能量函数E(f^t+1,V_t)，对应的水平集函数f^t+1的梯度下降流为：

$\begin{array}{c}\frac{\partial f^{t+1}}{\partial \tau }=\mu ·({\mathrm{\Delta f}}^{t+1}-div(\frac{▿f^{t+1}}{\left|▿f^{t+1}\right|}\left)\right)\\ +\kappa ·\left(\delta \right(f^{t+1})·div(\phi \_color·\frac{▿f^{t+1}}{\left|▿f^{t+1}\right|}\left)\right)\\ +\lambda ·\left(\right(1-in\_Silhouette\left(vox\right))·\delta (-f^{t+1})-in\_Silhouette(vox)·\delta (f^{t+1}\left)\right)\end{array}---\left(15\right)$

则水平集函数f^t+1的演化方程：

$f_{itt+1}^{t+1}(x,y,z)=f_{itt}^{t+1}(x,y,z)+\delta \tau ·\frac{\partial f_{itt}^{t+1}}{\partial t}---\left(16\right)$

其中，δτ表示水平集函数演化的时间步长；

将更新后的水平集函数作为已知条件，那么场景流的估计问题就相当于求取时空颜色相关性约束项中关于场景流的梯度下降流：

$\frac{\partial V_t^{vox}}{\partial t}=-\delta \left(f^{t+1}\right)·\frac{\partial \phi \_color}{\partial V_t^{vox}}·\left|▿f^{t+1}\right|---\left(17\right)$

场景流相应的估计过程为：

$V_t^{vox}=V_t^{vox}+\delta \tau ·\frac{\partial V_t^{vox}}{\partial \tau }---\left(18\right)$

当水平集函数停止演化时，其对应的零水平集即为时刻t+1的目标三维外形，此时的V_t即为所估计的场景流。

本发明的有益效果是：本发明多视环境中基于水平集函数演化的非刚性目标场景流估计，用水平集函数的演化模拟非刚性目标的运动，在时空颜色相关性约束、投影区域一致性约束下有效地驱动非刚性目标初始外形曲面向目标曲面演化，而水平集函数自身特点的约束，使的水平集函数演化后仍为水平集函数，避免了水平集函数重新初始化的问题，解决了现在估计非刚性目标的场景流时，由于环境因素，计算方法不当，往往会使得场景流不够连续、不平滑或存在空洞的问题。

具体实施方式

下面结合具体实施方式对本发明进行详细说明。

本发明多视环境中基于水平集函数演化的非刚性目标场景流估计，具体按照以下步骤实施：

步骤1：多视环境设置及预处理

将nCamNum个摄像机环绕观测场景放置，并使用张氏标定法获取摄像机的内、外参数；然后使用nCamNum个摄像机采集不包含目标的背景图像{Back₁,Back₂,...,Back_j,...,Back_nCamNum}；之后将动态非刚性目标O置于观测场景中，多个摄像机同步采集目标的运动视频，在时刻t、时刻t+1，同步采集到的视频帧分别为{I_1,t,I_2,t,...,I_j,t,...,I_nCamNum,t}、{I_1,t+1,I_2,t+1,...,I_j,t+1,...,I_nCamNum,t+1}，最后使用背景减法获取目标O的侧影轮廓，即目标图像，在时刻t和时刻t+1，分别使用多视图像帧减去背景图像，从而获取目标图像：{Silhouette_1,t,...,Silhouette_j,t,...,Silhouette_nCamNum,t}、{Silhouette_1,t+1,...,Silhouette_j,t+1,...,Silhouette_nCamNum,t+1}。

步骤2：获取静态目标三维外形

静态目标三维外形的获取过程即为动态目标场景流估计所需水平集函数的初始化过程。

将多视摄像机的公共观测区域离散化为一个三维体素空间Ω；

首先，使用{Silhouette_1,t,...,Silhouette_j,t,...,Silhouette_nCamNum,t}，基于体素着色法重建目标三维外形的体素模型Shape_Data_t，其中，

然后，重建目标三维外形，及初始化场景流估计所需的水平集函数：

Step1：初始化水平集函数φ(x)：在Ω内，靠近Ω边缘的区域构造一个足够大的封闭初始曲面Γ(0)，确保其能够包括Shape_Data的所有体素，构造初始水平集函数φ₀(x)：

其中c₀为一预先给定的正整数；

Step2：构造关于水平集函数φ(x)的能量函数E(φ)：

$E\left(\phi \right)=\alpha ·\underset{\Omega }{\int }g\left(x\right)\delta \left(\phi \right)\left|▿\phi \right|dx+\beta ·P\left(\phi \right)+\gamma ·A\left(\phi \right)---\left(2\right)$

其中，α、β和γ为比例系数，g(x)为在Ω中计算距离场为：g(x)＝dist(x,Shape_Data)，其中，x∈Ω，dist(x,Shape_Data)表示体素x到三维外形Shape_Data的最短距离，δ为Dirac函数；

P(φ)为惩罚项为：

$P\left(\phi \right)=\underset{\Omega }{\int }\rho \left(\right|▿\phi \left|\right)dx---\left(3\right)$

为：

$\rho \left(\right|▿\phi \left|\right)=\left(\begin{array}{cc}\frac{1}{{\left(2\pi \right)}^2}·(1-\cos (2\pi ·\left|▿\phi \right|\left)\right)& \left|▿\phi \right|\le 1\\ \frac{1}{2}·{\left(\right|▿\phi |-1)}^2& \left|▿\phi \right|>1\end{array}\right)---\left(4\right)$

A(φ)为加速项为：

$A\left(\phi \right)=\underset{\Omega }{\int }H(g-\eta )g\left(x\right)H(-\phi )dx---\left(5\right)$

H(x)是Heaviside函数，η是偏移程度，当g＜η时，加速项失效，从而可确保曲面不会穿过空洞；在曲面演化开始前，η赋值为以体素为单位的模型最大空洞的直径长度。

对于能量函数E(φ)，其相应的梯度下降流，即水平集函数的演化速度为：

Step3：令更新水平集函数：

φⁿ⁺¹(x)＝φⁿ(x)+Δτ·D(φⁿ(x))>

n表示迭代次数，Δτ为时间步长；

为了加快计算速度，本发明使用窄带法：提取零水平集φ(x,τ)＝0，将提取出的体素放入集合LevelSet中，对于每个体素i，i∈LevelSet，构造窄带其中为体素i的邻域，r为邻域的范围，通常指定其为以i为中心，2r+1为边长的立方体，即更新只在零水平集的邻域内进行。

Step4：判断是否小于阈值来判断是否收敛，若不收敛，则按照式(7)计算φⁿ⁺¹(x)的值，并转至Step4，若收敛，则结束；

令f^t(x)＝φⁿ⁺¹(x)，作为时刻t目标三维外形的水平集函数表示。

步骤3：构建时间序列上的关于水平集函数的能量泛函

将多视环境获取的非刚性目标运动视频提供的信息作为约束条件，构建关于水平集函数的能量泛函E(f^t+1,V_t)(f^t+1为水平集函数，其初始值为f^t，V_t表示时刻t目标三维外形的场景流)。在该能量函数的最小化过程中，水平集函数演化，其相应零水平集(f^t+1(vox)＝0时)向着运动后的目标三维外形逼近。

为了保证零水平集函数能够向着目标三维外形演化，本发明使用了如下三种约束条件：

①时空颜色相关性约束项

时空颜色相关性约束项用于驱动水平集函数向着目标三维外形演化，包括：时间颜色相关性和空间颜色相关性两部分。时间颜色相关性φ_temp(vox,V_t^vox)是指在时刻t目标三维外形上的某一体素以速度V_t^vox运动，在时刻t+1对应的目标三维外形上的体素为vox，其运动前后的颜色保持一致。空间颜色相关性φ_space(vox,t+1)是指时刻t+1目标三维外形上体素vox与其周围的体素具有颜色一致性。

时间颜色相关性φ_temp(vox,V_t^vox)为：

$\begin{array}{c}\phi \_temp(vox,V_t^{vox})\\ =\frac{1}{nVisNum}\underset{j=1}{\overset{nVisNum}{\Sigma }}\left(\frac{\partial photo\left({\mathrm{pro}}_j\right(vox),t+1)}{\partial vox}\right)·V_t^{vox}+\frac{\partial photo\left({\mathrm{pro}}_j\right(vox),t+1)}{\partial (t+1)}{)}^2\end{array}---\left(8\right)$

nVisNum表示对体素vox可见的视点个数，pro_j(vox)表示体素vox在视点j的二维投影，photo(X,t+1)为时刻t+1像素X对应的像素颜色值，φ_temp(vox,V_t^vox)反映了目标三维外形上的体素运动前后在多视图像中投影的颜色相关性：φ_temp(vox,V_t^vox)越小，颜色相关性越大；

空间颜色相关性φ_space(vox,t+1)为：

$\phi \_space(vox,t+1)=\frac{1}{nVisNum}\underset{j=1}{\overset{nVisNum}{\Sigma }}{\left(photo\right({\mathrm{pro}}_j\left(vox\right),t+1)-\overline{color})}^2---\left(9\right)$

其中：

$\overline{color}=\frac{1}{nVisNum}\underset{j=1}{\overset{nVisNum}{\Sigma }}photo\left({\mathrm{pro}}_j\right(vox),t+1)$

将时间颜色相关性和空间颜色相关性结合，得出时空颜色相关性约束项：

$\begin{array}{c}\phi \_color(vox,V_t^{vox})\\ =isColored(vox,V_t^{vox})·\phi \_temp(vox,V_t^{vox})+(1-isColored(vox,V_t^{vox}\left)\right)·\phi \_space(vox,t+1)\end{array}---\left(10\right)$

其中，isColored(vox,V_t^vox)表示时刻t体素vox-V_t^vox是否被判断为目标三维外形上的体素且已被着色，其取值方式为：

此时，Ω中所有体素的时空颜色相关性约束对应的能量函数为：

$E_{color}(f^{t+1},V_t)=\underset{vox\in \Omega }{\int \int \int }\phi \_color(vox,V_t^{vox})·\left|▿H\right|dxdydz---\left(11\right)$

②投影区域一致性约束项

投影区域一致性约束项用于避免水平集函数的过度演化。

目标三维外形在多视二维图像中的投影都会落入到该目标的侧影轮廓中：因此当表示目标三维外形的零水平集上的体素在所有视点的投影均落入到目标在各个视点的侧影轮廓中时，水平集函数的演化速度为正，即零水平集向着非目标区域运动；反之，水平集演化速度为负，即零水平集向着目标区域运动。

投影区域一致性约束项对应的能量函数为：

$E_{region}\left(f^{t+1}\right)=\underset{vox\in \Omega }{\int \int \int }(in\_Silhouette(vox)·H(f^{t+1})+(1-in\_Silhouette\left(vox\right))·H(-f^{t+1}\left)\right)dxdydz---\left(12\right)$

其中in_Silhouette(vox)表示体素vox在多视图像中的投影是否都落入目标的侧影轮廓中：

③水平集函数规整项

水平集函数规整项用于保持水平集函数在演化过程中的稳定性：保证水平集函数每次演化后仍为水平集函数，避免了水平集函数演化后的重新初始化问题。该约束项是由水平集函数的定义来确定的，即：由于水平集函数f为距离函数，而距离函数满足约束条件本发明也以该约束条件作为水平集函数的规整项。

水平集函数规整项对应的能量函数P(f^t+1)为：

$P\left(f^{t+1}\right)=\underset{\Omega }{\int }\rho \left(\right|▿f^{t+1}\left|\right)dx---\left(13\right)$

在上述三种约束条件下，构造时刻t+1关于水平集函数f^t+1和场景流V_t的能量泛函E(f^t+1,V_t)，该能量泛函由三部分组成：关于时空颜色相关性约束项的能量函数E_color(f^t+1,V_t)、关于投影区域一致性约束项的能量函数E_region(f^t+1)，水平集函数规整项的能量函数P(f^t+1)，即：

E(f^t+1,V_t)＝κ·E_color(f^t+1,V_t)+λ·E_region(f^t+1)+μ·P(f^t+1)>

其中，μ＞0为常量参数，用于控制f偏离符号距离函数后的惩罚效应，κ、λ表示权重系数。

步骤4：获取动态非刚性目标三维外形及其场景流

优化能量函数E(f^t+1,V_t)，获得相应的水平集函数f^t+1及对应场景流V_t。本发明基于梯度下降流法获取当能量函数的最小化时对应的水平集函数。f对应的零水平集即为运动后非刚性目标的三维外形，而此时就相当于在能量函数E(f^t+1,V_t)中，已知f^t+1，求取场景流V_t，具体为：

最小化能量函数E(f^t+1,V_t)，对应的水平集函数f^t+1的梯度下降流为：

$\begin{array}{c}\frac{\partial f^{t+1}}{\partial \tau }=\mu ·({\mathrm{\Delta f}}^{t+1}-div(\frac{▿f^{t+1}}{\left|▿f^{t+1}\right|}\left)\right)\\ +\kappa ·\left(\delta \right(f^{t+1})·div(\phi \_color·\frac{▿f^{t+1}}{\left|▿f^{t+1}\right|}\left)\right)\\ +\lambda ·\left(\right(1-in\_Silhouette\left(vox\right))·\delta (-f^{t+1})-in\_Silhouette(vox)·\delta (f^{t+1}\left)\right)\end{array}---\left(15\right)$

则水平集函数f^t+1的演化方程：

$f_{itt+1}^{t+1}(x,y,z)=f_{itt}^{t+1}(x,y,z)+\delta \tau ·\frac{\partial f_{itt}^{t+1}}{\partial t}---\left(16\right)$

其中，δτ表示水平集函数演化的时间步长；

将更新后的水平集函数作为已知条件，那么场景流的估计问题就相当于求取时空颜色相关性约束项中关于场景流的梯度下降流：

$\frac{\partial V_t^{vox}}{\partial t}=-\delta \left(f^{t+1}\right)·\frac{\partial \phi \_color}{\partial V_t^{vox}}·\left|▿f^{t+1}\right|---\left(17\right)$

场景流相应的估计过程为：

$V_t^{vox}=V_t^{vox}+\delta \tau ·\frac{\partial V_t^{vox}}{\partial \tau }---\left(18\right)$

当水平集函数停止演化时，其对应的零水平集即为时刻t+1的目标三维外形，此时的V_t即为所估计的场景流。