详解LK光流法(含金字塔多层光流),反向光流法(附代码)
LK光流法可用来跟踪特征点的位置。
比如在img1中的特征点,由于相机或物体的运动,在img2中来到了不同的位置。后面会称img1为Template(T),img2为I。
光流法有个假设:
灰度不变假设:同一个空间点的像素灰度值,在各图像中是不变的,也就是说T中特征点处的灰度,到了I中仍然是一样的灰度。
这就要求光照恒定,物体反射恒定,是个很强的假设。
现在要估计的是运动偏移量[dx, dy],也就是光流。仅用一个点无法解,一般会取一个窗口内的像素,考虑它们具有相同的运动。
用最小二乘法来解像素的运动如下:
这里的W表示对图像 I 进行像素变换,把它变到T所在的时间,即运动前的图像,看跟T有多大误差。这里假设对窗口内的像素进行变换。
p是运动位移量(dx, dy)。
这是个非线性优化问题,因为像素是跟坐标不是线性关系。
这时可以假设p已经知道了(假设就是运动前的点坐标,或者给定一个值),在它的基础上不断加上增量进行调整。
于是变成了优化如下问题,把它叫做error:
每次估计出增量Δp\Delta pΔp
更新p,
(4)(5)两步反复迭代直到达到收敛条件,收敛条件可以是
Δp\Delta pΔp 小于一个阈值,或者(4)中的cost比前一次大(理论上cost是逐渐减小的)
上面就是最小二乘法求光流的大概步骤,具体如何求Δp\Delta pΔp,下面是高斯牛顿法解Δp\Delta pΔp,及迭代出光流(dx, dy)的步骤:
把上面的(4)式,也就是cost函数,对I(W(x;p+Δp))I(W(x;p+\Delta p))I(W(x;p+Δp))进行一阶泰勒展开,得到:
这里的▽I\bigtriangledown I▽I指在图像 I,也就是img2中对特征点处求x,y方向上的梯度,然后通过W变换回T的坐标中。
∂W/∂p\partial W/ \partial p∂W/∂p指W对p求偏导。
举个例子吧,假如W如下
那么[∂Wx/∂p1∂Wx/∂p2∂Wy∂p1∂Wy/∂p2]=[1001]\left[ \begin{matrix} \partial W_x/ \partial p1 & \partial W_x/ \partial p2 \\ \partial W_y \partial p1 & \partial W_y/ \partial p2 \\ \end{matrix} \right] = \left[ \begin{matrix} 1 &0 \\ 0 & 1 \\ \end{matrix} \right] [∂Wx/∂p1∂Wy∂p1∂Wx/∂p2∂Wy/∂p2]=[1001]
上面的(p1, p2)就是我们要求的光流(dx, dy), 这个∂W/∂p\partial W/ \partial p∂W/∂p就是雅可比矩阵JJJ上面(6)式对Δp\Delta pΔp求偏导,得到:
先解释一下Δp\Delta pΔp和上面p的关系,上面的p是我们要求的光流(dx, dy),而直接求比较困难,现在假设已经知道初始的(dx, dy), 比如(0, 0),要通过每次求(dx, dy)的增量(Δdx,Δdy)\Delta dx, \Delta dy)Δdx,Δdy), 也就是Δp\Delta pΔp,来不断修正(dx, dy), 直到收敛。
另上面(9)式偏导为0,可求得:
这就是最小二乘法的解,
其实可以想到最小二乘法解增量问题的形式一般是
HΔx=bH \Delta x = bHΔx=b, 其中
H=JJT,b=−JeH = JJ^{T}, b = -JeH=JJT,b=−Je, eee即为error函数。
同样地,上面的H=JJTH = JJ^{T}H=JJT,
Δp\Delta pΔp一旦求出,就可以不断迭代,得到最终的光流(dx, dy), 也就是p了。
下面通过一个例子说明如何追踪特征点的光流。
先在img1,也就是T里提取特征点kp1
Mat img1 = imread("../imgs/LK1.png", 0); //T
Mat img2 = imread("../imgs/LK2.png", 0); //I//key points
vector<KeyPoint> kp1;
FAST(img1, kp1, 40); //可以用其他特征点
对每个特征点kp1[i], 设初始(dx, dy) = (0, 0)
auto kp = kp1[i];
double dx = 0, dy = 0; //需要估计
要用到的几个矩阵
Eigen::Matrix2d H = Eigen::Matrix2d::Zero(); //Hessian
Eigen::Vector2d b = Eigen::Vector2d::Zero(); //bias
Eigen::Vector2d J; //Jacobian
前面说了,不能只计算一个点,要计算一个小窗口内的点,并假设它们的运动是一样的。
代码中出现了逆向光流法,这个后面解释。
下面的代码中求出了H, b, 而Δp=H−1b\Delta p = H^{-1}bΔp=H−1b
//计算cost和jacobian
for(int x = -half_patch_size; x < half_patch_size; x++) {for(int y = -half_patch_size; y < half_patch_size; y++) {double error = GetPixelValue(img1, kp.pt.x + x, kp.pt.y + y) -GetPixelValue(img2, kp.pt.x + x + dx, kp.pt.y + y + dy);if(inverse == false) {//error分别对dx, dy求导,因为是离散的,所以用dx+1,dx-1的位移差,dy同样//x, y是窗口内的偏移量//img1中没有dx, dy的变量,所以只有img2J = -1.0 * Eigen::Vector2d(0.5 * (GetPixelValue(img2, kp.pt.x + x + dx + 1, kp.pt.y + y + dy) -GetPixelValue(img2, kp.pt.x + x + dx - 1, kp.pt.y + y + dy)),0.5 * (GetPixelValue(img2, kp.pt.x + x + dx, kp.pt.y + y + dy + 1) -GetPixelValue(img2, kp.pt.x + x + dx, kp.pt.y + y + dy - 1)));}else if (iter == 0) {//反向光流法,只需计算第一次的H,后面H固定J = -1.0 * Eigen::Vector2d(0.5 * (GetPixelValue(img1, kp.pt.x + x + 1, kp.pt.y + y) -GetPixelValue(img1, kp.pt.x + x - 1, kp.pt.y + y)),0.5 * (GetPixelValue(img1, kp.pt.x + x, kp.pt.y + y + 1) -GetPixelValue(img1, kp.pt.x + x, kp.pt.y + y - 1)));}//计算H, b and set costb += -error * J;cost += error * error;//正向光流法每次迭代都要更新Hif(inverse == false || iter == 0) {H += J * J.transpose();}}
Δp=H−1b\Delta p = H^{-1}bΔp=H−1b
Eigen::Vector2d update = H.ldlt().solve(b); //避免求矩阵的逆//最小二乘收敛
if(iter > 0 && cost > lastCost) {break;
}
迭代更新(dx, dy)
dx += update[0];
dy += update[1];
lastCost = cost;
最后(dx, dy)就是我们要求的光流,根据img1中的keypoint kp1, 可以追踪到img2中的keypoint坐标为
kp2[i].pt = kp.pt + Point2f(dx, dy);
以上是单层光流,下面说说金字塔的多层光流
假设左边最下面一层(3号)是原图像img1 (T), 那么往上(2号,1号)就是对img1的缩放,右边最下面一层是img2, 往上是对img2的缩放。
double pyramid_scale = 0.5;
//create pyramids
vector<Mat> pyr1, pyr2; //image pyramids
for(int i = 0; i < pyramids; i++) {if(i == 0) {pyr1.push_back(img1);pyr2.push_back(img2);} else {Mat img1_pyr, img2_pyr;//cv::Size(列,行)resize(pyr1[i - 1], img1_pyr,Size(pyr1[i-1].cols * pyramid_scale, pyr1[i-1].rows * pyramid_scale));resize(pyr2[i - 1], img2_pyr,Size(pyr2[i-1].cols * pyramid_scale, pyr2[i-1].rows * pyramid_scale));pyr1.push_back(img1_pyr);pyr2.push_back(img2_pyr);}
}
计算光流时,从最上面一层(1号)计算。
左边1号缩放的img1为T,右边1号缩放的img2为I,以img1的keypoint kp1计算光流,得到kp2。
到计算下面两层的时候,以上一层光流的结果作为下一层的初始假设,而不是像最上层那样假设(dx, dy) = (0, 0)。
double dx = 0, dy = 0; //需要估计
if(has_initial) {dx = kp2[i].pt.x - kp.pt.x; //x位移dy = kp2[i].pt.y - kp.pt.y; //y位移
}
这样做有什么好处呢?假设移动了20个像素,直接求光流可能由于变化太大而陷入局部极值;但是缩放之后可能就只移动了5个像素,这就是一种coarse to fine的思想。
现在来说逆向光流法,
前面正向光流每迭代一次都要计算一次H,计算量很大。考虑有没有一种方法可以只计算一次H,然后后面都用这个H。
逆向光流法是前面正向光流的逆向,也就是交换一下方向,现在是从I 反向回到T,也就是从运动之后的 I 变换回运动前的T。
引用一个图方便理解
error函数如下,只需要交换T 和 I
计算得到:
其中
可以看到由于T 是运动前的img,并没有移动p=(dx, dy),因此H和p无关,每次迭代计算Δp\Delta pΔp时H是恒定的,只需要在第一次迭代时计算一次即可。对应了上面代码的如下部分:
}else if (iter == 0) {//反向光流法,只需计算第一次的H,后面H固定//H只和T,也就是img1有关J = -1.0 * Eigen::Vector2d(0.5 * (GetPixelValue(img1, kp.pt.x + x + 1, kp.pt.y + y) -GetPixelValue(img1, kp.pt.x + x - 1, kp.pt.y + y)),0.5 * (GetPixelValue(img1, kp.pt.x + x, kp.pt.y + y + 1) -GetPixelValue(img1, kp.pt.x + x, kp.pt.y + y - 1)));}
每次迭代中H不需再更新
if(inverse == false || iter == 0) {H += J * J.transpose();
}
上面就是正向光流,反向光流,如果想了解更加详细的资料,请参考paper
完整版代码请参考链接
详解LK光流法(含金字塔多层光流),反向光流法(附代码)相关推荐
- MCMC详解2——MCMC采样、M-H采样、Gibbs采样(附代码)
MCMC是一种随机采样方法,用来处理一些复杂运算的近似求解.在HMM.LDA等模型中都有重要应用. 上一篇 MCMC详解1--蒙特卡洛方法 目录 1,马尔科夫链模型 1.1 马尔科夫链转移矩阵的性质 ...
- vue template html属性,详解template标签用法(含vue中的用法总结)
一.html5中的template标签 html中的template标签中的内容在页面中不会显示.但是在后台查看页面DOM结构存在template标签.这是因为template标签天生不可见,它设置了 ...
- 详解python实现FP-TREE进行关联规则挖掘(带有FP树显示功能)附源代码下载(3)
详解python实现FP-TREE进行关联规则挖掘(带有FP树显示功能)附源代码下载(3) 上一节简单讲了下FP树的生成,在这一节我将描述FP树的挖掘过程. 首先我们回顾一下要挖掘的特征项及样本空间: ...
- python中fp是什么意思_详解python实现FP-TREE进行关联规则挖掘(带有FP树显示功能)附源代码下载(3)...
详解python实现FP-TREE进行关联规则挖掘(带有FP树显示功能)附源代码下载(3) 上一节简单讲了下FP树的生成,在这一节我将描述FP树的挖掘过程. 首先我们回顾一下要挖掘的特征项及样本空间: ...
- python编程入门与案例详解pdf-这些年我读过的技术经典图书(附电子版下载地址)...
C技术资料 1.<> 作者: 谭浩强 这是我推荐的第一本书, 也是我接触的第一本书, 为什么把它放在第一位, 因为我觉得这本书对我的影响很大, 感觉这本书的最大特点是: 内容很全面, 内容 ...
- 操作指令详解_爱码小士丨 APP稳定性测试(附视频详解)
在实际的测试过程中,主要是对系统的功能来进行测试,用于校验功能的正确性 还需要考虑到系统在未修改的状态下,是否能够稳定运行,即崩溃.闪退.重启.系统异常等等等地情况 在APP中,稳定性测试一般是交由M ...
- 【Spring AOP】@Aspect结合案例详解(二): @Pointcut使用@within和within(已附源码)
文章目录 前言 @within 完善打印日志案例 @within深入说明 within 匹配指定类 匹配指定包(package) 源码下载 总结 前言 在微服务流行的当下,在使用Spring Clou ...
- simulink仿真实例详解_三菱FX 5U PLC模块硬件精品实例,附接线图
今天说说三菱FX5U 模块硬件的接线实例,主要有以下几个方面:电源AC.DC接线.输入输出接线.模拟量接线.不同原理有不同的接线方式,现在给大家仔细讲解分享! AC电源接线例 漏型输入[-公共端]时的 ...
- ❤️十大排序算法详解❤️——可能是你看过最全的,完整版代码
文章目录 前言 交集排序 冒泡 简单 快速排序 插入排序 直接插入排序 希尔排序 选择排序 简单选择排序 堆排序 归并排序 二路 多路 非比较类 计数排序 桶排序 基数排序 最后 前言 兄弟们,应上篇 ...
最新文章
- Jetty 基本使用样例
- Maven 打包的3中场景
- Java基础 Day04(个人复习整理)
- 逆向so_安卓逆向 | 分析调试与so调用实战
- CodeForces 362B 	Petya and Staircases
- Windows家庭版远程服务
- Laravel系列教程一:安装及环境配置
- [转贴]记那对住在我隔壁储藏室的大学刚毕业的小夫妻
- 第一章 WEB应用程序开发流程
- matlab控制电动机调速,控制电机调速及matlab仿真.doc
- perl对日志进行压缩备份小程序
- 验证IBIS模型正确并转换为DML模型图文演示
- unity3d:百度语音在线语音转文字,文字转语音,跨平台
- C语言实现windows窗口滑动条,四、Windows子窗口控件的滚动条类别—窗口子类别化(Window Subc...
- matplotlib animation动画保存(save函数)详解
- 目录即服务如何帮助安全、高效地管理WiFi用户?
- [WUSTCTF2020]佛说:只能四天
- Scanpy Umap 3D
- 太空射击第16课: 道具(Part 2)
- git reset 的三种模式的使用场景
热门文章
- 日程日历安排组件使用(angular)
- 张一鸣讲述字节跳动7年:诞生于民宅,大力出奇迹
- 《PHP从入门到精通》学习笔记之一
- c语言与指针——(二)指针变量的定义与赋值
- ViewPager设置setPageTransformer后RecyclerView垂直滑动问题
- OpenAPI报错集合
- python学习实验报告(第九周)
- Qt 使用C++特性“引用” - 获得槽函数的返回值
- 看到一个经典的魔兽版评论,玩过的可以进来看看,超级搞笑!!
- C++ 并行编程(thread)