Farneback 光流算法详解与 calcOpticalFlowFarneback 源码分析
来源:互联网 发布:linux mmap使用场景 编辑:程序博客网 时间:2024/05/22 01:30
光流基础概念
真实的三维空间中,描述物体运动状态的物理概念是运动场。三维空间中的每一个点,经过某段时间的运动之后会到达一个新的位置,而这个位移过程可以用运动场来描述。
而在计算机视觉的空间中,计算机所接收到的信号往往是二维图片信息。由于缺少了一个维度的信息,所以其不再适用以运动场描述。光流场(optical flow)就是用于描述三维空间中的运动物体表现到二维图像中,所反映出的像素点的运动向量场。
- 当描述部分像素时,称为:稀疏光流
- 当描述全部像素时,称为:稠密光流
光流法是利用图像序列中的像素在时间域上的变化、相邻帧之间的相关性来找到的上一帧跟当前帧间存在的对应关系,计算出相邻帧之间物体的运动信息的一种方法。光流法理解的关键点有:
- 核心问题:同一个空间中的点,在下一帧即将出现的位置
- 重要假设:光流的变化(向量场)几乎是光滑
- 角点处的光流能够通过角点的邻域完全确定下来,因此角点处的运动信息最为可靠;其次是边界的信息
光流法有着各种各样的分支,本文介绍的则是一种被广泛使用的经典稠密光流算法:Farneback 光流算法
Farneback 光流算法
- 作者:Gunnar Farneback
- 参考资料
- Two-Frame Motion Estimation Based on Polynomial Expansion,提纲阅读
- Polynomial Expansion for Orientation and Motion Estimation,作者博士毕业论文,深读(以下简称:博士论文)
- 源码:https://searchcode.com/file/30099587/opencv_source/src/cv/cvoptflowgf.cpp
OpenCV 主函数源码:
void calcOpticalFlowFarneback( const Mat& prev0, const Mat& next0, Mat& flow0, double pyr_scale, int levels, int winsize, int iterations, int poly_n, double poly_sigma, int flags ){ …… for( k = 0, scale = 1; k < levels; k++ ) { scale *= pyr_scale; if( prev0.cols*scale < min_size || prev0.rows*scale < min_size ) break; } levels = k; for( k = levels; k >= 0; k-- ) { for( i = 0, scale = 1; i < k; i++ ) scale *= pyr_scale; …… 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 ); } prevFlow = flow; }}
源代码中为了解决孔径问题,主函数中引入了图像金字塔。
孔径问题(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
OpenCV 子函数 FarnebackPolyExp:
static void FarnebackPolyExp( const Mat& src, Mat& dst, int n, double sigma )
理论基础
图像建模
将图像视为二维信号的函数(输出图像是灰度图像),因变量是二维坐标位置
其中,
因此系数化之后,以上公式等号右侧可以写为:
求解空间转换
如果将原有(笛卡尔坐标系)图像的二维信号空间,转换到以
Farneback 算法对于每帧图像中的每个像素点周围设定一个邻域
在一个邻域内灰度值的
此处博士论文中有举例说明,非常便于理解,详见博士论文 3.4 小节
权重分配
利用最小二乘法求解时,并非是邻域内每个像素点样本误差都对中心点有着同样的影响力,函数中利用二维高斯分布将影响力赋予权重。
在一个邻域内二维高斯分布的
对偶 (dual) 转换
为了“进一步加快求解”得到系数矩阵
而通过对偶转换之后,计算系数向量
函数输出
子函数 FarnebackPolyExp 输出得到的是单张图像中每个像素点的系数向量
Separable Normalized Convolution
博士论文 4.4 小节中提出的 Separable Normalized Convolution 计算方法提出将卷积操作由一维的直接计算拆分成两个维度的分别计算,可以降低计算复杂度,拆分的依据源自:
源码解读
- 【基于邻域】产生二维高斯分布的基础是一维高斯分布,一维高斯分布存储于数组 g 中,并且进行了求和后归一化:
double s = 0.;for( x = -n; x <= n; x++ ){ g[x] = (float)std::exp(-x*x/(2*sigma*sigma)); s += g[x];}s = 1./s;
- 【基于邻域】基于 Separable Normalized Convolution 方法,分别求解
(1,x,x2) 上的权重分布为:
for( x = -n; x <= n; x++ ){ g[x] = (float)(g[x]*s); xg[x] = (float)(x*g[x]); xxg[x] = (float)(x*x*g[x]);}
- 【基于邻域】求解对偶转换矩阵
G ,注意此时的二维高斯分布权重已经被按照列向量拉伸成为了一个(2n+1)2×1 的向量。而根据一维高斯分布数组 g 生成二维高斯分布权重也很简单,嵌套两层循环即可。循环时可以利用矩阵的对称性减少计算量。并且因为邻域中心点为(0,0) 点,所以循环求和时G(0,1)+=g[y]*g[x]*x
这类计算在x
和-x
上相互抵消,结果必然为 0 无需计算:
Mat_<double> G = Mat_<double>::zeros(6, 6);for( y = -n; y <= n; y++ ) for( x = -n; x <= n; x++ ) { G(0,0) += g[y]*g[x]; G(1,1) += g[y]*g[x]*x*x; G(3,3) += g[y]*g[x]*x*x*x*x; G(5,5) += g[y]*g[x]*x*x*y*y; }//G[0][0] = 1.;G(2,2) = G(0,3) = G(0,4) = G(3,0) = G(4,0) = G(1,1);G(4,4) = G(3,3);G(3,4) = G(4,3) = G(5,5);// invG:// [ x e e ]// [ y ]// [ y ]// [ e z ]// [ e z ]// [ u ]Mat_<double> invG = G.inv(DECOMP_CHOLESKY);double ig11 = invG(1,1), ig03 = invG(0,3), ig33 = invG(3,3), ig55 = invG(5,5);
注意:实际计算中
- 分别在 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); }}
OpenCV 子函数 FarnebackUpdateMatrices
static void FarnebackUpdateMatrices( const Mat& _R0, const Mat& _R1, const Mat& _flow, Mat& _M, int _y0, int _y1 )
子函数 FarnebackUpdateMatrices 中需要的变量参数包括:
1. _R0
:输入前一帧图像(编号:0)
2. _R1
:输入后一帧图像(编号:1)
3. _flow
:已知的前一帧图像光流场(用于 Blur 子函数的迭代,以及主函数中图像金字塔的不同层数间迭代)
4. _M
:储存中间变量,不断更新
5. _y0
:图像求解光流的起始行
6. _y1
:图像求解光流的终止行
理论基础
利用 FarnebackPolyExp 函数分别得到了视频前后帧中每个像素点的系数向量之后,则需要利用系数向量计算得到光流场。这个计算过程在博士论文 7.6 小节中有描述。
首先,每个像素点都有着初始位移(最开始设置为全 0 变量),将上一帧的初始位移增加到第一帧图像上的像素点位置
1. 将初始位移四舍五入得到整数值(博士论文中采用)
2. 利用二次线性插值(OpenCV 中的默认做法)
由此得到用于计算的中间变量
如果涉及到尺度变换,例如图像金字塔技术(子函数 FarnebackUpdateMatrices 中再次使用 multi-scale 的思路,但需要和主函数中的图像金字塔区分开来)。还需要另外与尺度缩放矩阵
源码解读
- 二次插值得到新一帧图像位置中的系数向量:
if( (unsigned)x1 < (unsigned)(width-1) && (unsigned)y1 < (unsigned)(height-1) ){ float a00 = (1.f-fx)*(1.f-fy), a01 = fx*(1.f-fy), a10 = (1.f-fx)*fy, a11 = fx*fy; // 2D interpolation r2 = a00*ptr[0] + a01*ptr[5] + a10*ptr[step1] + a11*ptr[step1+5]; r3 = a00*ptr[1] + a01*ptr[6] + a10*ptr[step1+1] + a11*ptr[step1+6]; r4 = a00*ptr[2] + a01*ptr[7] + a10*ptr[step1+2] + a11*ptr[step1+7]; r5 = a00*ptr[3] + a01*ptr[8] + a10*ptr[step1+3] + a11*ptr[step1+8]; r6 = a00*ptr[4] + a01*ptr[9] + a10*ptr[step1+4] + a11*ptr[step1+9]; r4 = (R0[x*5+2] + r4)*0.5f; r5 = (R0[x*5+3] + r5)*0.5f; r6 = (R0[x*5+4] + r6)*0.25f;}else{ r2 = r3 = 0.f; r4 = R0[x*5+2]; r5 = R0[x*5+3]; r6 = R0[x*5+4]*0.5f;}r2 = (R0[x*5] - r2)*0.5f;r3 = (R0[x*5+1] - r3)*0.5f;r2 += r4*dy + r6*dx;r3 += r6*dy + r5*dx;
- 处理不同尺度上的缩放,这里借鉴的是 multi-scale 思路,详见博士论文的 7.7 小节,目的是为了提高本算法的鲁棒性:
if( (unsigned)(x - BORDER) >= (unsigned)(width - BORDER*2) || (unsigned)(y - BORDER) >= (unsigned)(height - BORDER*2)){ float scale = (x < BORDER ? border[x] : 1.f)* (x >= width - BORDER ? border[width - x - 1] : 1.f)* (y < BORDER ? border[y] : 1.f)* (y >= height - BORDER ? border[height - y - 1] : 1.f); r2 *= scale; r3 *= scale; r4 *= scale; r5 *= scale; r6 *= scale;}
- 计算中间变量
G 与h ,存储到矩阵M
中:
M[x*5] = r4*r4 + r6*r6; // G(1,1)M[x*5+1] = (r4 + r5)*r6; // G(1,2)=G(2,1)M[x*5+2] = r5*r5 + r6*r6; // G(2,2)M[x*5+3] = r4*r2 + r6*r3; // h(1)M[x*5+4] = r6*r2 + r5*r3; // h(2)
OpenCV 子函数 FarnebackUpdateFlow_GaussianBlur、FarnebackUpdateFlow_Blur
static voidFarnebackUpdateFlow_GaussianBlur( const Mat& _R0, const Mat& _R1, Mat& _flow, Mat& _M, int block_size, bool update_matrices );static voidFarnebackUpdateFlow_Blur( const Mat& _R0, const Mat& _R1, Mat& _flow, Mat& _M, int block_size, bool update_matrices );
理论基础
根据光流的基本假设:光流的变化(向量场)几乎是光滑的。因此利用中间变量
源码解读
- 根据中间变量的元素,计算得到光流场存储于 flow 数组。flow 数组也是主函数最后的输出量:
double idet = 1./(g11*g22 - g12*g12 + 1e-3);flow[x*2] = (float)((g11*h2-g12*h1)*idet);flow[x*2+1] = (float)((g22*h1-g12*h2)*idet);
- 因为模糊化操作可能会执行多次,需要调用子函数 FarnebackUpdateMatrices 更新中间变量矩阵 M:
y1 = y == height - 1 ? height : y - block_size;if( update_matrices && (y1 == height || y1 >= y0 + min_update_stripe) ){ FarnebackUpdateMatrices( _R0, _R1, _flow, _M, y0, y1 ); y0 = y1;}
函数使用举例
#include "opencv2/highgui.hpp"#include "opencv2/imgproc.hpp"#include "opencv2/objdetect/objdetect.hpp"#include "opencv2/video/tracking.hpp"#include <vector>#include <stdio.h>#include <iostream>using namespace cv;using namespace std;int main(){ VideoCapture cap(0); Mat flow, frame; UMat flowUmat, prevgray; for(;;) { bool Is = cap.grab(); if(Is == false){ cout << "Video Capture Fail" << endl; break; } else{ Mat img; Mat original; cap.retrieve(img, CV_CAP_OPENNI_BGR_IMAGE); resize(img, img, Size(640, 480)); img.copyTo(original); cvtColor(img, img, COLOR_BGR2GRAY); if(prevgray.empty() == false){ calcOpticalFlowFarneback(prevgray, img, flowUmat, 0.4, 1, 12, 2, 8, 1.2, 0); flowUmat.copyTo(flow); for(int y=0; y<original.rows; y+= 5){ for(int x=0; x<original.cols; x+= 5){ const Point2f flowatxy = flow.at<Point2f>(y, x)*10; line(original, Point(x,y), Point(cvRound(x+flowatxy.x), cvRound(y+flowatxy.y)), Scalar(255, 0, 0)); circle(original, Point(x,y), 1, Scalar(0,0,0), -1); } } namedWindow("prew", WINDOW_AUTOSIZE); imshow("prew", original); img.copyTo(prevgray); } else{ img.copyTo(prevgray); } int key1 = waitKey(20); } }}
- Farneback 光流算法详解与 calcOpticalFlowFarneback 源码分析
- 【光流】Farneback\brox\SF\Flownet2光流算法比较
- Opencv Python版学习笔记(四)光流跟踪之Gunnar Farneback’s 算法
- Farneback 稠密光流--求两幅图像之间的光流--代码(本人略有修改)
- calcOpticalFlowFarneback
- calcOpticalFlowFarneback
- OpenCV中CalcOpticalFlowFarneback()函数分析
- OpenCV中CalcOpticalFlowFarneback()函数分析
- 基于光流分析的运动目标快速检测与跟踪融合算法
- 计算两幅图片的farneback 稠密光流,并将结果图显示出来的程序
- 【算法分析】Lucas–Kanade光流算法
- 【算法分析】Lucas–Kanade光流算法
- 【算法分析】Lucas–Kanade光流算法
- sift算法与源码分析
- LK光流算法
- 光流算法
- 光流算法
- LK光流算法
- 深入理解java异常处理机制
- bash shell(bash) 快捷键
- JAVA实现一个简单的代数运算语言编译器(二)--词法分析准备
- 学习日记
- express无法在页面中无法获取本地图片的解决方法
- Farneback 光流算法详解与 calcOpticalFlowFarneback 源码分析
- jdbc连接Oracle数据库
- Python self参数 & 函数详解
- <<摸着石头过河>>摘录四
- 北京圣思园Java教学视频全集迅雷下载
- Programming Ability Test 1003.我要通过
- 水<并查集,但是WA好久>
- 蓝桥杯之高精度算法
- Git Rebase教程: 用Git Rebase让时光倒流