真实的三维空间中,描述物体运动状态的物理概念是运动场。三维空间中的每一个点,经过某段时间的运动之后会到达一个新的位置,而这个位移过程可以用运动场来描述。
而在计算机视觉的空间中,计算机所接收到的信号往往是二维图片信息。由于缺少了一个维度的信息,所以其不再适用以运动场描述。光流场(optical flow)就是用于描述三维空间中的运动物体表现到二维图像中,所反映出的像素点的运动向量场。
当描述部分像素时,称为:稀疏光流当描述全部像素时,称为:稠密光流光流法是利用图像序列中的像素在时间域上的变化、相邻帧之间的相关性来找到的上一帧跟当前帧间存在的对应关系,计算出相邻帧之间物体的运动信息的一种方法。光流法理解的关键点有:
核心问题:同一个空间中的点,在下一帧即将出现的位置重要假设:光流的变化(向量场)几乎是光滑角点处的光流能够通过角点的邻域完全确定下来,因此角点处的运动信息最为可靠;其次是边界的信息光流法有着各种各样的分支,本文介绍的则是一种被广泛使用的经典稠密光流算法:Farneback 光流算法
源代码中为了解决孔径问题,主函数中引入了图像金字塔。
孔径问题(Aperture Problem): - http://blog.csdn.net/hankai1024/article/details/23433157 - 形象理解:从小孔中观察一块移动的黑色幕布观察不到任何变化。但实际情况是幕布一直在移动中 - 解决方案:从不同尺度(图像金字塔)上对图像进行观察,由高到低逐层利用上一层已求得信息来计算下一层信息
主函数 calcOpticalFlowFarneback 中需要的变量参数包括: 1. _prev0:输入前一帧图像 2. _next0:输入后一帧图像 3. _flow0:输出的光流 4. pyr_scale:金字塔上下两层之间的尺度关系 5. levels:金字塔层数 6. winsize:均值窗口大小,越大越能 denoise 并且能够检测快速移动目标,但会引起模糊运动区域 7. iterations:迭代次数 8. poly_n:像素邻域范围大小,一般为 5、7 等 9. poly_sigma:高斯标准差,一般为 1~1.5(函数处理中需要高斯分布权重) 10. flags:计算方法,包括 OPTFLOW_USE_INITIAL_FLOW 和 OPTFLOW_FARNEBACK_GAUSSIAN
实际的 Farneback 算法在每一层金字塔上的处理过程为:
Mat R[2], I, M; for( i = 0; i < 2; i++ ) { img[i]->convertTo(fimg, CV_32F); GaussianBlur(fimg, fimg, Size(smooth_sz, smooth_sz), sigma, sigma); resize( fimg, I, Size(width, height), CV_INTER_LINEAR ); FarnebackPolyExp( I, R[i], poly_n, poly_sigma ); } FarnebackUpdateMatrices( R[0], R[1], flow, M, 0, flow.rows ); for( i = 0; i < iterations; i++ ) { if( flags & OPTFLOW_FARNEBACK_GAUSSIAN ) FarnebackUpdateFlow_GaussianBlur( R[0], R[1], flow, M, winsize, i < iterations - 1 ); else FarnebackUpdateFlow_Blur( R[0], R[1], flow, M, winsize, i < iterations - 1 ); }输入的图像默认为灰度图片,默认只有亮度值
图像输入与输出时进行的高斯模糊操作(Gaussian Blur)操作都是使得光流场结果平滑,从而满足假设:光流的变化几乎是光滑的
所以需要关注的子函数有 4 个: 1. FarnebackPolyExp 2. FarnebackUpdateMatrices 3. FarnebackUpdateFlow_GaussianBlur 4. FarnebackUpdateFlow_Blur
将图像视为二维信号的函数(输出图像是灰度图像),因变量是二维坐标位置 x=(xy)T ,并且利用二次多项式对于图像进行近似建模的话,会得到:
f(x)∼xTAx+bTx+c 其中, A 是一个2 × 2的矩阵, b 是一个2 × 1的矩阵因此系数化之后,以上公式等号右侧可以写为:
(xy)T(r4r6/2r6/2r5)(xy)+(r2r3)T(xy)+r1=r1+r2x+r3y+r4x2+r5y2+r6xy如果将原有(笛卡尔坐标系)图像的二维信号空间,转换到以 (1,x,y,x2,y2,xy) 作为基函数的空间,则表示图像需要一个六维向量作为系数,代入不同像素点的位置 x,y 求出不同像素点的灰度值。
Farneback 算法对于每帧图像中的每个像素点周围设定一个邻域 (2n+1)×(2n+1) ,利用邻域内的共 (2n+1)2 个像素点作为最小二乘法的样本点,拟合得到中心像素点的六维系数。因此对于图像中的每个像素点,都能得到一个六维向量。
在一个邻域内灰度值的 (2n+1)×(2n+1) 矩阵中,将矩阵按列优先次序拆分组合成 (2n+1)2×1 的向量 f ,同时已知 (1,x,y,x2,y2,xy) 作为基函数的转换矩阵 B 维度为 (2n+1)2×6 (也可以视为 6 个列向量 bi 共同组成的矩阵),邻域内共有的系数向量 r 维度为 6×1 ,则有:
f=B×r=(b1b2b3b4b5b6)×r此处博士论文中有举例说明,非常便于理解,详见博士论文 3.4 小节
利用最小二乘法求解时,并非是邻域内每个像素点样本误差都对中心点有着同样的影响力,函数中利用二维高斯分布将影响力赋予权重。 在一个邻域内二维高斯分布的 (2n+1)×(2n+1) 矩阵中,将矩阵按列优先次序拆分组合成 (2n+1)2×1 的向量 a 。因此原本的基函数的转换矩阵 B 将变为:
B=(a⋅b1a⋅b2a⋅b3a⋅b4a⋅b5a⋅b6)为了“进一步加快求解”得到系数矩阵 r ,博士论文 4.3 小节中提出使用对偶的方式再次转换基函数矩阵 B ,此时的对偶转换矩阵为 G ,经过转换后的基函数矩阵的列向量为 bi~ 。博士论文中 G 的计算方式为:
G=⎛⎝⎜⎜⎜(a⋅b1,b1)⋮(a⋅b6,b1)⋯⋱⋯(a⋅b1,b6)⋮(a⋅b6,b6)⎞⎠⎟⎟⎟ 而通过对偶转换之后,计算系数向量 r 的方式就简单了很多,其中 ⋆ 表示两个信号的 互相关(实质上类似于两个函数的卷积 https://zh.wikipedia.org/wiki/互相关 )过程: r(x)=G−1⎛⎝⎜⎜⎜((a⋅b1)⋆f)(x)⋮((a⋅b6)⋆f)(x)⎞⎠⎟⎟⎟子函数 FarnebackPolyExp 输出得到的是单张图像中每个像素点的系数向量 r (不包括常数项系数 r1 ,因为之后的计算光流过程中没有用到)
博士论文 4.4 小节中提出的 Separable Normalized Convolution 计算方法提出将卷积操作由一维的直接计算拆分成两个维度的分别计算,可以降低计算复杂度,拆分的依据源自:
注意:实际计算中 G 的特殊结构会使得 G−1 的特殊结构,所以只需要保存逆矩阵中几个特殊位置的元素即可。证明过程见博士论文的附录 A
分别在 vertical 和 horizontal 两个方向进行卷积计算,卷积结果分别存在 row 数组和 b1~b6 中,最终的系数向量存在于 drow 数组中(不需要保存常量系数 r1,因为不涉及后面子函数的计算过程)。计算时注意 x 和 -x 在计算时导致的正负性差异: for( y = 0; y < height; y++ ) { float g0 = g[0], g1, g2; float *srow0 = (float*)(src.data + src.step*y), *srow1 = 0; float *drow = (float*)(dst.data + dst.step*y); // vertical part of convolution for( x = 0; x < width; x++ ) { row[x*3] = srow0[x]*g0; row[x*3+1] = row[x*3+2] = 0.f; } for( k = 1; k <= n; k++ ) { g0 = g[k]; g1 = xg[k]; g2 = xxg[k]; srow0 = (float*)(src.data + src.step*std::max(y-k,0)); srow1 = (float*)(src.data + src.step*std::min(y+k,height-1)); for( x = 0; x < width; x++ ) { float p = srow0[x] + srow1[x]; float t0 = row[x*3] + g0*p; float t1 = row[x*3+1] + g1*(srow1[x] - srow0[x]); // (x) and (-x) float t2 = row[x*3+2] + g2*p; row[x*3] = t0; row[x*3+1] = t1; row[x*3+2] = t2; } } // horizontal part of convolution // row padding: left & right for( x = 0; x < n*3; x++ ) { row[-1-x] = row[2-x]; row[width*3+x] = row[width*3+x-3]; } for( x = 0; x < width; x++ ) { g0 = g[0]; // r1 ~ 1, r2 ~ x, r3 ~ y, r4 ~ x^2, r5 ~ y^2, r6 ~ xy double b1 = row[x*3]*g0, b2 = 0, b3 = row[x*3+1]*g0, b4 = 0, b5 = row[x*3+2]*g0, b6 = 0; for( k = 1; k <= n; k++ ) { double tg = row[(x+k)*3] + row[(x-k)*3]; g0 = g[k]; b1 += tg*g0; // 1 * 1 b4 += tg*xxg[k]; // 1 * x^2 b2 += (row[(x+k)*3] - row[(x-k)*3])*xg[k]; // 1 * x b3 += (row[(x+k)*3+1] + row[(x-k)*3+1])*g0; // y * 1 b6 += (row[(x+k)*3+1] - row[(x-k)*3+1])*xg[k]; // y * x b5 += (row[(x+k)*3+2] + row[(x-k)*3+2])*g0; // y^2 * 1 } // do not store r1 drow[x*5+1] = (float)(b2*ig11); drow[x*5] = (float)(b3*ig11); drow[x*5+3] = (float)(b1*ig03 + b4*ig33); drow[x*5+2] = (float)(b1*ig03 + b5*ig33); drow[x*5+4] = (float)(b6*ig55); } }子函数 FarnebackUpdateMatrices 中需要的变量参数包括: 1. _R0:输入前一帧图像(编号:0) 2. _R1:输入后一帧图像(编号:1) 3. _flow:已知的前一帧图像光流场(用于 Blur 子函数的迭代,以及主函数中图像金字塔的不同层数间迭代) 4. _M:储存中间变量,不断更新 5. _y0:图像求解光流的起始行 6. _y1:图像求解光流的终止行
利用 FarnebackPolyExp 函数分别得到了视频前后帧中每个像素点的系数向量之后,则需要利用系数向量计算得到光流场。这个计算过程在博士论文 7.6 小节中有描述。
首先,每个像素点都有着初始位移(最开始设置为全 0 变量),将上一帧的初始位移增加到第一帧图像上的像素点位置 x 上,得到此像素点在下一帧图像上的大致位置 x̃ :
x̃ =x+d̃ (x) x̃ 可能是浮点位置,则无法估计此位置的系数向量值。对待此问题的方法有两种: 1. 将初始位移四舍五入得到整数值(博士论文中采用) 2. 利用二次线性插值(OpenCV 中的默认做法)由此得到用于计算的中间变量 A(x),Δb(x) :
A(x)=A1(x)+A2(x̃ )2Δb(x)=−12(b2(x̃ )−b1(x))+A(x)d̃ (x) 如果涉及到尺度变换,例如图像金字塔技术(子函数 FarnebackUpdateMatrices 中再次使用 multi-scale 的思路,但需要和主函数中的图像金字塔区分开来)。还需要另外与尺度缩放矩阵 S(x) 有关。总之,子函数 FarnebackUpdateMatrices 的目的是为了得到中间变量 G 与 h : G(x)=S(x)TA(x)TA(x)S(x)h(x)=S(x)TA(x)TΔb(x)根据光流的基本假设:光流的变化(向量场)几乎是光滑的。因此利用中间变量 G 与 h 求解光流场 dout 前,需要进行一次局部模糊化处理(主函数中 winsize 输入变量控制),可选均值模糊 (FarnebackUpdateFlow_Blur)、高斯模糊 (FarnebackUpdateFlow_GaussianBlur)。对于模糊后的中间变量,可以直接求解光流场:
dout(x)=Gavg(x)−1havg(x)