---

欢迎访问个人网络日志🌹🌹知行空间🌹🌹

---

光流法介绍

光流描述了像素在图像中的运动，就像彗星☄划过天空中流动图像。同一个像素，随着时间的流逝，会在图像中运动，光流法就是追踪它的运动过程。

光流法根据追踪的像素数又可以分成稀疏光流法和稠密光流法。

• 稀疏光流法：计算部分像素的运动，稀疏法以Lucas-Kanade 光流为代表，可以用来目标追踪中跟踪特征点的位置。
• 稠密光流法：计算所有像素的运动，稠密光流法以Horn-Schunck光流为代表。

在Lucas-Kanade光流中，将相机的图像看成是随时间变化的，图像$I$在$t$时刻位置为$(x,y)$处的像素，它的灰度值可以写成：

$$I(x,y,t)$$

通过这种方式将图像看成了关于位置和时间的函数。

考虑固定的空间点，在世界坐标系中的其位置是固定不变的，在$t$时刻，其在图像中的像素坐标为$(x,y)$。由于相机在运动，该空间点在$t+1$时刻在图像中的像素坐标将发生变动，如何估计在$t+1$时刻同个空间点的像素坐标呢？这正是光流法要解决的问题。

光流法的基本假设：同一个空间点的像素灰度值，在各个图像中的是固定不变的。

上述假设使用公式描述就是说，

$t$时刻在图像中位置$(x,y)$处的像素 $I(x,y,t)$

在$t+dt$时刻运动到了图像的$(x+dx,y+dy)$处，

基于灰度假设有以下关系成立:

$$I(x,y,t)=I(x+dx,y+dy,t+dt)$$

灰度假设是一个很强的假设，实际中很可能不成立，由于物体的材质/相机成想的角度/光照条件发生变化的时候，同一个空间点的像素灰度值很有可能发生变化，因此光流法的结果不一定可靠。

在此假设成立的前提下，来看下如何计算像素的运动。

对上式右侧进行泰勒展开：

$$I(x+dx,y+dy,t+dt)\approx I(x,y,t)+\frac{\partial I}{\partial x}dx+\frac{\partial I}{\partial y}dy+\frac{\partial I}{\partial t}dt$$

因为假设了灰度不变，即$I(x+dx,y+dy,t+dt)=I(x,y,t)$,因此：

$$\frac{\partial I}{\partial x}dx+\frac{\partial I}{\partial y}dy+\frac{\partial I}{\partial t}dt=0$$

对上式两边同时除以$dt$得:

$$\frac{\partial I}{\partial x}\frac{dx}{dt}+\frac{\partial I}{\partial y}\frac{dy}{dt}=-\frac{\partial I}{\partial t}$$

$\frac{dx}{dt}$为像素在x轴上的速度，记为$\mu$

$\frac{dy}{dt}$为像素在y轴上的速度，记为$v$

$\frac{\partial I}{\partial x}$为图像在该点$x$方向上的梯度，记为$I_x$

$\frac{\partial I}{\partial y}$为图像在该点$y$方向上的梯度，记为$I_y$

图像灰度值对时间的变化量，记为$I_t$（？_?不是已经假设了图像灰度值不变的吗？_?）

写成矩阵形式有：

$$\left[\begin{array}{cc}{I}_x& I_y\end{array}\right]\left[\begin{array}{c}\mu \\ v\end{array}\right]=-I_t$$

上式中，

$I_x/I_y$是x/y方向的图像梯度，可以从图像中直接求得。

$I_t$是像素灰度值的变化量，也可从图像中直接求得。

$\mu /v$是像素在x,y方向的运动速度，有了速度，就可以求距离了，因此$\mu /v$就是我们期望计算的量。

而上式是二元一次方程，单纯的通过上式还无法求出$\mu /v$。

在LK光流中，引入了新的假设作为约束，即：某一个窗口内的像素具有相同的运动。

考虑一个大小为$w\times w$的窗口，其包含$w^2$个像素。因为假设该窗口内像素具有相同的运动，因此可以得到$w^2$个方程,

$${\left[\begin{array}{cc}{I}_x& I_y\end{array}\right]}_k\left[\begin{array}{c}\mu \\ v\end{array}\right]=-I_{tk},k=1,2,\mathrm{3..}{w}^2$$

记：
$$A=\left[\begin{array}{c}[I_x,I_y{]}_1\\ ...\\ [I_x,I_y{]}_k\end{array}\right],b=\left[\begin{array}{c}{I}_{t1}\\ ...\\ I_{tk}\end{array}\right]$$

上面的方程就变成了，

$$A\left[\begin{array}{c}\mu \\ v\end{array}\right]=b$$

这个是关于$\mu /v$的超定线性方程，可以使用最小二乘法求解。

$${\left[\begin{array}{c}\mu \\ v\end{array}\right]}^{\ast }=-(A^{T}A{)}^{-1}{A}^{T}b$$

OpenCV中calcOpticalFlowPyrLK函数

该方法使用迭代Lucas-Kanade算法计算稀疏特征点的光流，用来做特征点跟踪，该方法使用了金字塔，因此具有一定的尺度不变性。

void cv::calcOpticalFlowPyrLK(	InputArray 	        prevImg,InputArray 	        nextImg,InputArray 	        prevPts,InputOutputArray 	nextPts,OutputArray 	    status,OutputArray 	    err,Size 	            winSize = Size(21, 21),int 	            maxLevel = 3,TermCriteria 	    criteria = TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 30, 0.01),int 	            flags = 0,double          	minEigThreshold = 1e-4 
);

• prevImg:上一帧图像
• nextImg:下一帧图像
• prevPts:上一帧图像中关键点
• nextPts:根据光流计算的上一帧关键点在当前帧中的位置
• status:关键点的跟踪状态，vector,1表示OK,0表示LOST
• err:每个特征点的跟踪误差vector<float>有len(status)==len(nextPts)==len(prePts)==len(err)
• winSize:在每层金字塔中，LK算法中用来求解计算像素运动而假设具有相同运动的窗口大小。
• maxLevel:金字塔的层数，层数多，尺度不变性能更好，运算时间更久
• criteria:迭代搜索算法的终止条件,默认值表示在指定的最大迭代次数criteria.maxCount(30)之后或当搜索窗口移动小于criteria.epsilon(0.01)时终止迭代
• flags:设置误差或者初始值参数,可选下面两个值：
  • OPTFLOW_USE_INITIAL_FLOW设置使用nextPts中的值作为迭代的初始值，如果不设置为OPTFLOW_USE_INITIAL_FLOW,初始状态就使用prevPts中的值，直接从prevPts复制到nextPts,OpenCV源码中对OPTFLOW_USE_INITIAL_FLOW的使用方式为:

    if( flags & OPTFLOW_USE_INITIAL_FLOW )nextPt = nextPts[ptidx]*(float)(1./(1 << level));
    elsenextPt = prevPt;
    

    • OPTFLOW_LK_GET_MIN_EIGENVALS,flags设置为这个值时使用光流运动方程2x2的正规矩阵，也即空间梯度矩阵的最小特征值作为误差项。如果不设置成OPTFLOW_LK_GET_MIN_EIGENVALS,将原始点和移动点周围像素的$L_1$距离除以窗口中的像素作为误差项。
• minEigThreshold:迭代LK算法会计算光流运动方程2x2的正规矩阵，也即空间梯度矩阵的最小特征值，然后再除以运动不变窗口中的像素总数作为一个误差评价标准，当其小于minEigThreshold时，说明这个点已经追踪不到了，会将其从追踪特征点中移除，避免其对应相素运动的计算，可提升性能。

calcOpticalFlowPyrLK通常和goodFeatureToTrack方法一起使用，先使用GFTTDetector提取特征点的位置，再使用calcOpticalFlowPyrLK追踪其在连续视频流中的位置，避免了特征描述子的计算和特征点的匹配，可以极大的提升追踪的性能。

#include <memory>
#include <vector>
#include <cstdlib>#include <opencv2/features2d.hpp>
#include <opencv2/opencv.hpp>class TestOpticalFlowLK {public:typedef std::shared_ptr<TestOpticalFlowLK> Ptr;TestOpticalFlowLK();~TestOpticalFlowLK() = default;void track(std::vector<cv::String> &filenames) const;private:cv::Ptr<cv::GFTTDetector> gftt_ptr_;};TestOpticalFlowLK::TestOpticalFlowLK()
{gftt_ptr_ = cv::GFTTDetector::create(500, 0.2, 50);
}void TestOpticalFlowLK::track(std::vector<cv::String> &filenames) const
{assert(filenames.size() > 1);std::vector<cv::KeyPoint> kps1;std::vector<cv::Point2f>pts1, pts2;std::vector<cv::Scalar> colors;cv::Mat last_img = cv::imread(filenames[0], 0), cur_img;cv::Mat mask(last_img.size(), CV_8UC1, 255);gftt_ptr_->detect(last_img, kps1, mask);for(auto &kp : kps1) {int r = (int)(255. * rand() / (RAND_MAX + 1.f));int g = (int)(255. * rand() / (RAND_MAX + 1.f));int b = (int)(255. * rand() / (RAND_MAX + 1.f));std::cout << "r:" << r << "g:" << g << "b:" << b << std::endl;colors.emplace_back(r, g, b);pts1.push_back(kp.pt);pts2.push_back(kp.pt);}std::vector<uchar> status;// cv::Mat err;std::vector<float> err;cv::cvtColor(mask, mask, cv::COLOR_GRAY2BGR);cv::Mat frame;for (auto &filename : filenames){std::cout << "filename: " << filename << std::endl;cur_img = cv::imread(filename, 0);cv::calcOpticalFlowPyrLK(last_img,cur_img,pts1,pts2,status,err,cv::Size(13, 13),3,cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 30, 0.01),cv::OPTFLOW_USE_INITIAL_FLOW);cv::cvtColor(cur_img, cur_img, cv::COLOR_GRAY2BGR);int cnt = 0;for(size_t i = 0; i < status.size(); i++) {std::cout << " " << err[i];if(!status[i]) continue;if(abs((pts1[i].x - pts2[i].x)) > 80 | abs((pts1[i].y - pts2[i].y)) > 80) continue;cv::line(mask, pts1[i], pts2[i], colors[i], 2);cv::circle(cur_img, pts2[i], 10, colors[i], 1);pts1[i].x = pts2[i].x;pts1[i].y = pts2[i].y;cnt += 1;}std::cout << std::endl;cv::addWeighted(mask, 0.5, cur_img, 0.5, -65, frame);cv::imshow("frame", frame);cv::waitKey(0);cv::cvtColor(cur_img, cur_img, cv::COLOR_BGR2GRAY);last_img = cur_img;cv::imwrite("frame.png", frame);}
}

上述代码的运行效果为：

完整演示代码和数据可以在以下仓库中找到：

https://gitee.com/lx_r/basic_cplusplus_examples

补充

因为假设了像素灰度值不变，还可以将其看成非线性优化问题，求解如下方程:

$$\begin{array}{c}\underset{\Delta x,\Delta y}{\min }||I_1(x,y)-I_2(x+\Delta x,y+\Delta y)|{|}_2^2\end{array}$$

从零开始实现LK光流在视觉SLAM十四讲中高博已经实现过了，更读详细信息可以参考slambook2仓库。

https://github.com/gaoxiang12/slambook2

---

欢迎访问个人网络日志🌹🌹知行空间🌹🌹

---

reference

致谢：理论部分来自于《视觉SLAM十四讲(第二版)》P201-210.

• 1.https://zhuanlan.zhihu.com/p/384651830
• 2.https://docs.opencv.org/4.6.0/dc/d6b/group__video__track.html#ga473e4b886d0bcc6b65831eb88ed93323
• 3.Pyramidal implementation of the lucas kanade feature tracker