GMCM2017-前景目标提取

最近独自一人参加完成了GMCM2017-D题，实际上是做几个问题，本人在机器视觉上虽然饶有兴趣，但时毕竟没有经验，所以做出来的效果也是差强人意，不过第三问做出来后真的是很开心。

可以简短地将问题概括如下：

对静态背景进行前景提取
对动态视频进行前景提取
对晃动视频进行前景提取
对多视角的视频进行有效的前景提取
前景提取的作用（群体性事件的预警？）

对静态背景进行前景提取

对于静态的背景，我们可以采用高斯混合模型（以下推理部分主要来自[1]）利用视频每一帧的信息拟合出一个背景，然后利用背景差分的办法提取出前景目标。记每一帧的图像中某一个点的像素值为序列X={X1,X2,...,Xt}X=\{X_1,X_2,...,X_t\}，其中XjX_j是相互独立的随机统计过程。

由于固定场景中，图中像素点相对稳定，假设每个像素点均服从正态分布，一般地，我们可以用KK个正态分布去描述一个像素点，KK个状态中一部分表示像素点为背景时候的状态，一部分表示像素点为前景时候的状态。为了减少计算量，一般KK值取3。

由于某个像素点的灰度值用随机变量XtX_t表示，其概率密度函数可以表示为：

f(Xt=x)=∑i=1kωi,tη(x,μi,t,Σi,t),

\begin{equation} f\left( {{X_t} = x} \right) = \sum\limits_{i = 1}^k {{\omega _{i,t}}} \eta \left( {x,{\mu _{i,t}},{\Sigma _{i,t}}} \right), \end{equation}

其中η(x,μi,t,Σi,t)\eta \left( {x,{\mu _{i,t}},{\Sigma _{i,t}}} \right)是tt时刻的第ii个正态分布，其均值为μi,t{{\mu _{i,t}}}，方差为Σi,t{{\Sigma _{i,t}}}，为第ii个正态分布在tt时刻的权重，并有∑i=1kωi,t=1\sum\limits_{i = 1}^k {{\omega _{i,t}}} = 1，其中：

η(x,μi,t,Σi,t)=1(2π)n/2∣∣Σi,t∣∣1/2⋅exp(−12(xi−μi,t)TΣ−1i,t(xi−μi,t)),

\begin{equation} \eta \left( {x,{\mu _{i,t}},{\Sigma _{i,t}}} \right) = \frac{1}{{{{\left( {2\pi } \right)}^{n/2}}{{\left| {{\Sigma _{i,t}}} \right|}^{1/2}}}} \cdot \exp \left( { - \frac{1}{2}{{\left( {{x_i} - {\mu _{i,t}}} \right)}^T}\Sigma _{i,t}^{ - 1}\left( {{x_i} - {\mu _{i,t}}} \right)} \right), \end{equation}
上式中 nn表示XtX_t的维数，由于本文算法在提取前景目标的时候采用的是灰度图像，所以我们可以取 n=1n=1。首先初始化模型的参数，设 μ0\mu_{0}表示第一帧所有像素点的平均灰度值， σ20\sigma_0^2为其对应的平均均方差，即：

μ0=1N∑t=0NXt,σ20=1N∑t=0N(Xt−μ0)2,

\begin{equation} {\mu _0} = \frac{1}{N}\sum\limits_{t = 0}^N {{X_t}} ,\begin{array}{*{20}{c}} {}\end{array} = \frac{1}{N}\sum\limits_{t = 0}^N {{{\left( {{X_t} - {\mu _0}} \right)}^2}} , \end{equation}

上式中NN表示视频的帧数，初始化ωi=1/k\omega_i=1/k。接下来，为了区分前景和背景图像，我们需要定义一个匹配函数χ(Xi,t)\chi \left( {{X_{i,t}}} \right)：

χ(Xi,t)={1,∣∣Xi,t−μi,t∣∣≤Dσi,t−10,∣∣Xi,t−μi,t∣∣>Dσi,t−1,

\begin{equation} \chi \left( {{X_{i,t}}} \right) = \left\{ \begin{array}{l} 1,\left| {{X_{i,t}} - {\mu _{i,t}}} \right| \le D{\sigma _{i,t-1}}\\ 0,\left| {{X_{i,t}} - {\mu _{i,t}}} \right| > D{\sigma _{i,t-1}} \end{array} \right., \end{equation}
若 χ(Xi,t)=1\chi \left( {{X_{i,t}}} \right)=1，我们称当前像素值 Xi,tX_{i,t}与其模板匹配，否则称为不匹配。上式中 DD为一个阈值参数，一般取2.5，μi,t\mu_{i,t}表示 Xi,tX_{i,t}的第 ii个正态分布的均值，σi,t−1\sigma _{i,t-1}类似地定义。如果 Xi,tX_{i,t}与 KK个正态分布中存在匹配的话，说明Xi,tX_{i,t}满足一个或多个正态分布的描述，匹配的分布我们按下式对分布进行更新，若不匹配的分布保存原值：

μi,t=(1−ρ)μi,t−1+ρXi,t,σ2i,t=(1−ρ)σ2i,t−1+ρ(Xi,t−μi,t−1)T(Xi,t−μi,t−1),ωi,t,k=(1−α)ωi,t−1,k+αMi,t,k,

\begin{array}{l} {\mu _{i,t}} = \left( {1 - \rho } \right){\mu _{i,t - 1}} + \rho {X_{i,t}},\\ \sigma _{i,t}^2 = \left( {1 - \rho } \right)\sigma _{i,t - 1}^2 + \rho {\left( {{X_{i,t}} - {\mu _{i,t - 1}}} \right)^T}\left( {{X_{i,t}} - {\mu _{i,t - 1}}} \right),\\ {\omega _{i,t,k}} = \left( {1 - \alpha } \right){\omega _{i,t - 1,k}} + \alpha {M_{i,t,k}}, \end{array}

其中α\alpha表示更新率，α\alpha越大更新速率越快，但是为了减少引入更多的背景噪声，一般选取较小的值，ρ\rho表示参数的更新率，定义为ρ=α/ωi,t\rho=\alpha/\omega_{i,t}，对于匹配的正态分布Mi,t,kM_{i,t,k}取为1，否则取为0，。如果像素点Xi,tX_{i,t}和所有的正态分布均不匹配，说明Xi,tX_{i,t}很可能是运动目标出现的点，或者运动目标移动过后显现的背景点，此时我们对KK个正态分布中权重最小的那个分布的均值用Xi,tX_{i,t}替代，另外初始化一个较大的标准差σ\sigma，其余正态分布保证不变，然后权重按上式权重更新公式进行更新。(取Mi,t,k=0M_{i,t,k}=0)

考虑一个视频图像，由于我们考察的背景是相对静止的。所以，对于每一个背景点，其方差是比较小的，根据上述更新公式可以知道，用于描述背景点的高斯分布的权重不断增大，考虑这两个特点，我们可以建立如下前景提取办法，首先对权重进行归一化处理：

ωi,t,k=ωi,t∑i=1kωi,t,i=1,2,...,k

\begin{equation} {\omega _{i,t,k}} = \frac{{{\omega _{i,t}}}}{{\sum\limits_{i = 1}^k {{\omega _{i,t}}} }},i = 1,2,...,k \end{equation}

接下来，根据ωi,t,k/σi,t,k\omega_{i,t,k}/\sigma_{i,t,k}，从大到小对每个像素的KKK个正态分布进行排序。注意到，权重越大，方差越小的话，比值\omega_{i,t,k}/\sigma_{i,t,k}ωi,t,k/σi,t,k\omega_{i,t,k}/\sigma_{i,t,k}越大，即成为背景点的可能越大。我们设置一个区分点sss，令排序前sss个正态分布作为背景模型，剩下的作为前景模型，

\begin{equation} s = \mathop {\arg \min }\limits_k \left\{ {\sum\limits_{i = 1}^k {{\omega _{i,t}} > T} } \right\}, \end{equation}

s=argmink{∑i=1kωi,t>T},

\begin{equation} s = \mathop {\arg \min }\limits_k \left\{ {\sum\limits_{i = 1}^k {{\omega _{i,t}} > T} } \right\}, \end{equation}

其中TTT为一个自定义的阈值，一般经验取值为0.85，最后，对于目标视频，我们对每个像素点的X_{i,t}Xi,tX_{i,t}与前sss个描述背景的正态分布进行匹配，如果存在一个匹配，那么令该点为背景点，否则令其为前景点。

我们利用本节模型对pedestrian视频进行处理，得到如下结果：

结果一般，差强人意。上图分别为视频在第18帧到第21帧的分离结果，其中图(a)是分离结果，图(b)是视频原图。可以看出距离摄像机较近的部分分离效果良好，较远的目标分离效果欠佳，总体而言可以较好地提取前景目标。

对动态视频进行前景提取

对于动态背景的视频，我们不妨假设背景的变化幅度是较小的。这样我们可以把它看作是一个微小的扰动（黔驴技穷），我们可以利用图像的形态学处理和中值滤波对背景进行处理，结合本节的高斯混合模型对waterSurface视频处理可以得到如下结果：

上图中(a)为直接用本节的高斯混合模型得到的分离结果，图(b)为进行形态学处理和中值滤波之后的结果，可以看出图(b)的结果基本消除了背景的动态扰动，具有较好的结果。

晃动视频的前景目标提取

我们不妨假设晃动的幅度不大，此时我们可以将这样的幅度变换近似地看成是仿射变换。如下图所示：

将视角统一后，我们便可以把晃动视频的情况转换为固定视角的视频，由此可以由上一节的高斯混合模型得出前景分离结果。下面我们将以目标视频的第一帧作为基准，将其他帧的视频全部矫正为第一帧的视角，注意到的是我们考虑的视频必须是短时间的，这样能保证，其他帧到第一帧的仿射变换总是存在的(更进一步地，实际上只要保证对应的像素点至少有两对即可)。

我们设基准图像上的像素为(x′,y′)T(x',y')^T，其它帧上对应的像素点为(x,y)T(x,y)^T，于是有如下仿射变换：

⎛⎝⎜x′y′1⎞⎠⎟=⎛⎝⎜a11a210a12a220b1b21⎞⎠⎟⎛⎝⎜xy1⎞⎠⎟,

\left( \begin{array}{l} x'\\ y'\\ 1 \end{array} \right) = \left( {\begin{array}{*{20}{c}}{{a_{11}}}}}}\\{{a_{21}}}}}}\\0&0&1\end{array}} \right)\left( \begin{array}{l} x\\ y\\ 1 \end{array} \right),

利用对视角变化等仿射变换具有强大稳定性的SIFT算法我们可以获得足够多的特征点，下面我们利用获得的足够多的特征点来拟合求解出仿射变换对应的变换矩阵。以变换矩阵的第一行为例，我们有a11x+a12yb1=x′a_{11}x+a_{12}yb_1=x'，若有NN组特征点，那么有：

a11x(j)+a12y(j)+b(j)1=x′(j),j=1,2,...,N

\begin{equation} a_{11}x^{(j)}+a_{12}y^{(j)}+b_1^{(j)}=x'^{(j)},j=1,2,...,N \end{equation}

上式可以简写为Xa=cXa=c，其中a=(a11,a12,b1)Ta=(a_{11},a_{12},b_1)^T，以及有X=(x(j),y(j))∈RN×2X=(x^{(j)},y^{(j)}) \in R^{N \times 2}，c=x′(j)c=x'^{(j)},由于这一般是一个超定方程，由最小二乘法的思想我们可以将问题转换为求解mina∥Xa−c∥22\mathop {\min }\limits_a \left\| {Xa - c} \right\|_2^2。

令E=∥Xa−c∥22E=\left\| {Xa - c} \right\|_2^2，展开得E=aTXTXa−2aTXTc+cTcE={a^T}{X^T}Xa - 2{a^T}{X^T}c + {c^T}c，下面对EE做关于aa的偏导有:

∇aE=2(XTXa−XTc)=0⇒a=(XTX)−1XTc

\begin{equation} \begin{array}{l} {\nabla _a}E = 2\left( {{X^T}Xa - {X^T}c} \right) = 0\\ \Rightarrow a = {\left( {{X^T}X} \right)^{ - 1}}{X^T}c \end{array} \end{equation}

分析可知，当a=(XTX)−1XTca = {\left( {{X^T}X} \right)^{ - 1}}{X^T}c的时候，EE可以取得最小值，注意到XTXX^TX可逆的概率为1，即使不可逆也可以重新获取一次特征点。类似地，我们可以求解完整的仿射变换矩阵aa。

获得上图后不能马上用高斯混合模型进行处理，这是因为可以看到校正后的视频(图(b))有大量的黑边，如果直接套用高斯混合模型，将会把黑边的变化造成的噪声也识别出来。我们必须找到一个合适的视角区域，由于我们可以认为处理的视频是一个较短的视频，基于这个假设，我们可以取所有校正后的视频的公共重叠部分。可以知道，由假设，这个公共重叠部分是足够大的。

获得重叠的视角区域后，我们对所有校正后的视频都进行裁剪，接着利用高斯混合模型可以获得前景目标。最后我们利用在校正时候获得的变换矩阵aa，重新构造一个仿射逆变换a−1a^{-1}，最终将前景目标的晃动复原，结果如下图所示：

多视角的前景目标提取

多个视角的视频能给我们更多的信息，我们可以根据这些信息对前景目标进行三维重建。如下图所示，我们令绿色区域为我们的三维空间中的观测区域，以三个摄像机为例。考虑观测空间中的一个体像素vv，若它和前景目标相交，则称为存在，否则称为不存在。用如下概率表示：

P{Ev(i)=1|I},

\begin{equation} P\left\{ {{E_v}\left( i \right) = 1|I} \right\}, \end{equation}
其中ii表示第ii个体像素，II是所有观测视频的即可，即I={I1,I2,...,IN}I=\{I_1,I_2,...,I_N\}，其中的IkI_k表示第kk个观测视频图像。

对P{Ev(i)=1|I}P\left\{ {{E_v}\left( i \right) = 1|I} \right\}，由贝叶斯公式有：

P{Ev(i)=1|I}=P{I|Ev(i)=1}P{Ev(i)=1}P{I}=P{Ev(i)=1}∏k=1NP{Ik}⋅∏k=1NP{Ik|Ev(i)=1}=P{Ev(i)=1}∏k=1NP{Ik}⋅∏k=1NP{Ev(i)=1|Ik}P{Ik}P{Ev(i)=1}=1P{Ev(i)=1}N−1⋅∏k=1NP{Ev(i)=1|Ik}

\begin{equation} \begin{aligned} P\left\{ {{E_v}\left( i \right) = 1|I} \right\}&=\frac{{P\left\{ {I|{E_v}\left( i \right) = 1} \right\}P\left\{ {{E_v}\left( i \right) = 1} \right\}}}{{P\left\{ I \right\}}}\\ &=\frac{{P\left\{ {{E_v}\left( i \right) = 1} \right\}}}{{\prod\limits_{k = 1}^N {P\left\{ {{I_k}} \right\}} }} \cdot \prod\limits_{k = 1}^N {P\left\{ {{I_k}|{E_v}\left( i \right) = 1} \right\}}\\ &=\frac{{P\left\{ {{E_v}\left( i \right) = 1} \right\}}}{{\prod\limits_{k = 1}^N {P\left\{ {{I_k}} \right\}} }} \cdot \prod\limits_{k = 1}^N {\frac{{P\left\{ {{E_v}\left( i \right) = 1|{I_k}} \right\}P\left\{ {{I_k}} \right\}}}{{P\left\{ {{E_v}\left( i \right) = 1} \right\}}}}\\ &=\frac{1}{{P{{\left\{ {{E_v}\left( i \right) = 1} \right\}}^{N - 1}}}} \cdot \prod\limits_{k = 1}^N {P\left\{ {{E_v}\left( i \right) = 1|{I_k}} \right\}} \end{aligned} \end{equation}

注意到1/P{Ev(i)=1}N−11/{P{{\left\{ {{E_v}\left( i \right) = 1} \right\}}^{N - 1}}}是空间本身前景目标的存在与否决定的，与观测无关，是一个常数，换句话说，我们有：

P{Ev(i)=1|I}∝∏k=1NP{Ev(i)=1|Ik},

\begin{equation} P\left\{ {{E_v}\left( i \right) = 1|I} \right\} \propto \prod\limits_{k = 1}^N {P\left\{ {{E_v}\left( i \right) = 1|{I_k}} \right\}} , \end{equation}

由此，如上图所示，假设我们知道三个摄像机的相对位置，我们就可以连结摄像机和体像素在每个IkI_k上做投影得到如图所示的红色区域，我们用红色区域中非0的值的像素数量与红色区域的像素数量的比值来近似表示P{Ev(i)=1|Ik}{P\left\{ {{E_v}\left( i \right) = 1|{I_k}} \right\}}。由此，我们可以遍历所有的像素体并计算P{Ev(i)=1|I}P\left\{ {{E_v}\left( i \right) = 1|I} \right\}，设置一个阈值TT，令：

χ(i)={1,P{Ev(i)=1|I}≥T0,P{Ev(i)=1|I}<T,

\begin{equation} \chi \left( i \right) = \left\{ \begin{array}{l} 1,P\left\{ {{E_v}\left( i \right) = 1|I} \right\} \ge T\\ 0,P\left\{ {{E_v}\left( i \right) = 1|I} \right\}

由此即可完成三维重建，得到一个信息更加完善的前景目标检测模型。

群体性事件的自动判别

由于是我一个人在比赛，实在没有时间，所以这部分说实话，纯水系列可以这么说，ε=(´ο｀*)))唉

监控视频对群体性事件的自动监控是非常有意义的一项研究，本文以人群逃散为例，来研究这个问题。我们假设人群逃散的时候，人群的形态、速度大致相同，那么我们可以考虑一种平均的情况，即只需要研究一个人逃散状态的变化特点即可。

考虑全局前景目标的重心变化，对观测视频利用第一节中的高斯混合模型进行前景提取后，由于只需要观测重心的变化，所以我们可以对全局做一个模板足够大的中值滤波，这样的处理可以得到更为准确的重心(去掉噪声点)，然后对tt时刻视频中所有前景目标的位置求得其重心(Gx,t,Gy,t)(G_{x,t},G_{y,t})。由于，两帧之间的时间间隔相同，所以重心的移动距离sts_t即可以表征人群的平均移动速度，用下式表示：

st=(Gx,t−Gx,t−1)2+(Gy,t−Gy,t−1)2−−−−−−−−−−−−−−−−−−−−−−−−−−−√,

\begin{equation} s_t=\sqrt{(G_{x,t}-G_{x,t-1})^2+(G_{y,t}-G_{y,t-1})^2}, \end{equation}

我们也可以设置一个阈值TT，当st>T,t∈[t1,t2]s_t>T,t\in [t_1,t_2]时表示[t1,t2][t_1,t_2]时间段可能发生人群逃散情况。下面考虑两种情况来验证模型，首先考虑的是逃散方向和摄像机在同一直线上的情况：

对于逃散方向和摄像机共线的情况，由于共线关系，随之时间的推移，前景点数量大大增加，而且重心前期的位置变化也会较为集中，从下图图(a)可以看到，前期时，重心在一个小区域内密集变化。目标奔跑时，重心位置显著变化，由图(b)可以得到在大约65帧时候开始奔跑，与原视频情况吻合。

对于逃散方向和摄像机垂直的情况，由于垂直关系，随之时间的推移，前景点数量变化，而且重心的位置变化较为分散，从下图图(a)可以看到，目标奔跑时，重心位置显著变化，由图(b)可以得到在大约37帧时候开始奔跑，与原视频情况吻合。

Reference

[1]刘静. 固定场景下的运动目标检测算法与实现[D]. 湖南大学, 2010.

[2]徐剑, 丁晓青, 王生进,等. 多视角多行人目标检测、定位与对应算法[J]. 清华大学学报(自然科学版), 2009(8):1139-1143.

[3]梁铭富, 房少梅, 黄中展,等. 基于突发事件影响力传播的双向人流疏散仿真[J]. 计算机应用, 2017, 37(5):1496-1502.

[4]塞利斯基艾海舟, 兴军亮. 计算机视觉：算法与应用[M]. 清华大学出版社, 2012.

[5]Mikolajczyk K, Schmid C. A performance evaluation of local descriptors[J]. IEEE Transactions on Pattern Analysis Machine Intelligence, 2005, 27(10):1615.

[6]Fleuret F, Berclaz J, Lengagne R, et al. Multicamera people tracking with a probabilistic occupancy map[J]. IEEE Transactions on Pattern Analysis Machine Intelligence, 2008, 30(2):267.

代码可以到github上找https://github.com/dedekinds/GMCM2017