b站视频： https://www.bilibili.com/video/BV1DE411a7uH?from=search&seid=6718989481068612540
光流（一）--综述概览
10626 
收藏 12光流除了提供远近外，还可以提供角度信息。与咱们的眼睛正对着的方向成90度方向运动的物体速度要比其他角度的快，当小到0度的时候，也就是物体朝着我们的方向直接撞过来，我们就是感受不到它的运动（光流）了，看起来好像是静止的。当它离我们越近，就越来越大（当然了，我们平时看到感觉还是有速度的，因为物体较大，它的边缘还是和我们人眼具有大于0的角度的）。
光流的概念是Gibson在1950年首先提出来的。它是空间运动物体在观察成像平面上的像素运动的瞬时速度，是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系，从而计算出相邻帧之间物体的运动信息的一种方法。一般而言，光流是由于场景中前景目标本身的移动、相机的运动，或者两者的共同运动所产生的。
当人的眼睛观察运动物体时，物体的景象在人眼的视网膜上形成一系列连续变化的图像，这一系列连续变化的信息不断“流过”视网膜（即图像平面），好像一种光的“流”，故称之为光流（optical flow）。光流表达了图像的变化，由于它包含了目标运动的信息，因此可被观察者用来确定目标的运动情况。
研究光流场的目的就是为了从图片序列中近似得到不能直接得到的运动场。运动场，其实就是物体在三维真实世界中的运动；光流场，是运动场在二维图像平面上（人的眼睛或者摄像头）的投影。
那通俗的讲就是通过一个图片序列，把每张图像中每个像素的运动速度和运动方向找出来就是光流场。那怎么找呢？咱们直观理解肯定是：第t帧的时候A点的位置是(x1, y1)，那么我们在第t+1帧的时候再找到A点，假如它的位置是(x2,y2)，那么我们就可以确定A点的运动了：(ux, vy) = (x2, y2) - (x1,y1)。
那怎么知道第t+1帧的时候A点的位置呢？ 这就存在很多的光流计算方法了。
1981年，Horn和Schunck创造性地将二维速度场与灰度相联系，引入光流约束方程，得到光流计算的基本算法。人们基于不同的理论基础提出各种光流计算方法，算法性能各有不同。Barron等人对多种光流计算技术进行了总结，按照理论基础与数学方法的区别把它们分成四种：基于梯度的方法、基于匹配的方法、基于能量的方法、基于相位的方法。近年来神经动力学方法也颇受学者重视。
1 基于梯度的方法
基于梯度的方法又称为微分法，它是利用时变图像灰度（或其滤波形式）的时空微分（即时空梯度函数）来计算像素的速度矢量。由于计算简单和较好的结果，该方法得到了广泛应用和研究。典型的代表是Horn-Schunck的光流计算方法，该方法在光流基本约束方程的基础上附加了全局平滑假设，从而计算出光流场。基于此思想，大量的改进算法不断提出。Nagel采用有条件的平滑约束，即通过加权矩阵的控制对梯度进行不同平滑处理；Black和Anandan针对多运动的估计问题，提出了分段平滑的方法。虽然很多基于梯度的光流估计方法取得了较好的光流估计，但由于在计算光流时涉及到可调参数的人工选取、可靠性评价因子的选择困难，以及预处理对光流计算结果的影响，在应用光流对目标进行实时检测与自动跟踪时仍存在很多问题。
2 基于匹配的方法
基于匹配的光流计算方法包括基于特征和区域的两种。基于特征的方法不断地对目标主要特征进行定位和跟踪，对目标大的运动和亮度变化具有鲁棒性（robustness）。存在的问题是光流通常很稀疏，而且特征提取和精确匹配也十分困难。基于区域的方法先对类似的区域进行定位，然后通过相似区域的位移计算光流。这种方法在视频编码中得到了广泛的应用。然而，它计算的光流仍不稠密。另外，这两种方法估计亚像素精度的光流也有困难，计算量很大。在考虑光流精度和稠密性时，基于匹配的方法适用。
3 基于能量的方法
基于能量的方法首先要对输入图像序列进行时空滤波处理，这是一种时间和空间整合。对于均匀的流场，要获得正确的速度估计，这种时空整合是非常必要的。然而，这样做会降低光流估计的空间和时间分辨率。尤其是当时空整合区域包含几个运动成分（如运动边缘）时，估计精度将会恶化。此外，基于能量的光流技术还存在高计算负荷的问题。此方法涉及大量的滤波器，目前这些滤波器是主要的计算消费。然而，可以预期，随着相应硬件的发展，在不久的将来，滤波将不再是一个严重的限制因素，所有这些技术都可以在帧速下加以实现。
4 基于相位的方法
Fleet和Jepson首次从概念上提出了相位信息用于光流计算的问题。因为速度是根据带通滤波器输出的相位特性确定的，所以称为相位方法。他们根据与带通速度调谐滤波器输出中的等相位轮廓相垂直的瞬时运动来定义分速度。带通滤波器按照尺度、速度和定向来分离输入信号。
基于相位的光流技术的综合性能是比较好的：速度估计比较精确且具有较高的空间分辨率，对图像序列的适用范围也比较宽。同时，这里仍有几个问题值得讨论：
(1)与基于能量的光流技术一样，基于相位的模型既有一定的生物合理性，又有较高的计算复杂性；
(2)尽管相位技术用两帧图像就可计算光流，但要获得足够的估计精度，就必须有一定的整合时间，这个延迟将会降低边缘处运动估计的时间分辨率；
(3)Fleet和Jespon的方法对输入图像序列中的时间混叠比较敏感。
5 神经动力学方法
计算机视觉研究的初衷就是为了模仿人类视觉系统的功能。然而人类理解与识别图像的能力与计算机形成了巨大的反差。视觉科学家们迫切期望借鉴人类处理图像的方法，以摆脱困境。对于光流计算来讲，如果说前面的基于能量或相位的模型有一定的生物合理性的话，那么近几年出现的利用神经网络建立的视觉运动感知的神经动力学模型则是对生物视觉系统功能与结构的更为直接的模拟。
Grossberg等人的视觉运动感知神经动力学模型描述了运动感知中视皮层简单细胞、复杂细胞、超复杂细胞以及视网膜双极细胞之间的相互作用，揭示了运动分割与组合、竞争与合作的神经整合机制。这个称为运动边界轮廓系统的神经网络解释了复杂运动图形上的局部模糊运动如何被积极地组织成一个整体一致的运动信号，给出了整体小孔问题的一个解。这个模型对于整体运动方向的判别非常有效，然而它却不能给出运动速度的大小。
Fay和Waxman模仿视网膜中的时空处理和大脑的视觉运动通路，基于并联动力学提出了一个多层神经网络，它涉及光适应、边缘增强和边缘速度提取等几个处理阶段。网络中，每个节点的动力学特性类似于具有可变电导的细胞膜，光适应利用神经元间的抑制作用来获取，空间对比度增强借助于一个修正的on-中心/off-周边反馈网络来实现，最后的速度估计由一个称为对传活化法的动力学方程来提取。这个神经网络在一并行机上实现了30帧/秒的帧速下的速度提取。遗憾的是它仅能提供运动边缘的法向速度估计，为了恢复整个模式的光流场，还必须用速度泛函方法将估计的法向流整合成一个致密的光流场。尽管用这些神经动力学模型来测量光流还很不成熟，然而这些方法及其结论为进一步研究打下了良好的基础，是将神经机制引入运动计算方面所做的极有意义的尝试。一并行机上实现了30帧/秒的帧速下的速度提取。遗憾的是它仅能提供运动边缘的法向速度估计，为了恢复整个模式的光流场，还必须用速度泛函方法将估计的法向流整合成一个致密的光流场。
目前，对光流的研究方兴未艾，新的计算方法还在不断涌现。这里对光流技术的发展趋势与方向提出以下几点看法：
(1)现有技术各有自己的优点与缺陷，方法之间相互结合，优势互补，建立光流计算的多阶段或分层模型，是光流技术发展的一个趋势;
(2)通过深入的研究发现，现有光流方法之间有许多共通之处。如微分法和匹配法的前提假设极为相似；某些基于能量的方法等效于区域匹配技术；而相位方法则将相位梯度用于法向速度的计算。这些现象并不是偶然的。Singh指出，现有各种光流估计方法基本上可以统一在一个框架之中，这个框架将光流信息分成两类：保持信息和邻域信息，光流场的恢复通过两种信息的提取和融合来实现。光流计算的统一框架的研究是这个领域的又一趋势；
(3)尽管光流计算的神经动力学方法还很不成熟，然而对它的研究却具有极其深远的意义。随着生物视觉研究的不断深入，神经方法无疑会不断完善，也许光流计算乃至计算机视觉的根本出路就在于神经机制的引入。神经网络方法是光流技术的一个发展方向。
其他的咱们先不说了，回归应用吧（呵呵，太高深了，自己说不下去了）。OpenCV中实现了不少的光流算法。
1）calcOpticalFlowPyrLK
通过金字塔Lucas-Kanade 光流方法计算某些点集的光流（稀疏光流）。理解的话，可以参考这篇论文：”Pyramidal Implementation of the Lucas Kanade Feature TrackerDescription of the algorithm”
2）calcOpticalFlowFarneback
用Gunnar Farneback 的算法计算稠密光流（即图像上所有像素点的光流都计算出来）。它的相关论文是："Two-Frame Motion Estimation Based on PolynomialExpansion"
3）CalcOpticalFlowBM
通过块匹配的方法来计算光流。
4）CalcOpticalFlowHS
用Horn-Schunck 的算法计算稠密光流。相关论文好像是这篇：”Determining Optical Flow”
5）calcOpticalFlowSF
这一个是2012年欧洲视觉会议的一篇文章的实现："SimpleFlow: A Non-iterative, Sublinear Optical FlowAlgorithm"，工程网站是：http://graphics.berkeley.edu/papers/Tao-SAN-2012-05/ 在OpenCV新版本中有引入。
稠密光流需要使用某种插值方法在比较容易跟踪的像素之间进行插值以解决那些运动不明确的像素，所以它的计算开销是相当大的。而对于稀疏光流来说，在他计算时需要在被跟踪之前指定一组点（容易跟踪的点，例如角点），因此在使用LK方法之前我们需要配合使用cvGoodFeatureToTrack()来寻找角点，然后利用金字塔LK光流算法，对运动进行跟踪。但个人感觉，对于少纹理的目标，例如人手，LK稀疏光流就比较容易跟丢。
至于他们的API的使用说明，我们直接参考OpenCV的官方手册就行：
IJCV2011有一篇文章，《A Database and Evaluation Methodology for Optical Flow》里面对主流的光流算法做了简要的介绍和对不同算法进行了评估。网址是：
http://vision.middlebury.edu/flow/
感觉这个文章在光流算法的解说上非常好，条例很清晰。想了解光流的，推荐看这篇文章。另外，需要提到的一个问题是，光流场是图片中每个像素都有一个x方向和y方向的位移，所以在上面那些光流计算结束后得到的光流flow是个和原来图像大小相等的双通道图像。那怎么可视化呢？这篇文章用的是Munsell颜色系统来显示。
关于孟塞尔颜色系统（MunsellColor System），可以看wikibaike。它是美国艺术家阿尔伯特孟塞尔（Albert H. Munsell，1858－1918）在1898年创制的颜色描述系统。
孟塞尔颜色系统的空间大致成一个圆柱形：
南北轴=明度（value），从全黑（1）到全白（10）。
经度=色相（hue）。把一周均分成五种主色调和五种中间色：红(R)、红黄(YR)、黄(Y)、黄绿(GY)、绿(G)、绿蓝(BG)、蓝(B)、蓝紫(PB)、紫(P)、紫红(RP)。相邻的两个位置之间再均分10份，共100份。
距轴的距离=色度（chroma），表示色调的纯度。其数值从中间（0）向外随着色调的纯度增加，没有理论上的上限（普通的颜色实际上限为10左右，反光、荧光等材料可高达30）。由于人眼对各种颜色的的敏感度不同，色度不一定与每个色调和明度组合相匹配。
具体颜色的标识形式为：色相+明度+色度。
在上面的那个评估的网站有这个从flow到color显示的Matlab和C++代码。但是感觉C++代码分几个文件，有点乱，然后我自己整理成两个函数了，并配合OpenCV的Mat格式。
下面的代码是用calcOpticalFlowFarneback来计算稠密光流并且用这个颜色系统来显示的。这个计算稠密光流的方法与其他几个相比还是比较快的，640x480的视频我的是200ms左右一帧，但其他的一般都需要一两秒以上。结果图中，不同颜色表示不同的运动方向，深浅就表示运动的快慢了。
void calcOpticalFlowFarneback(InputArray prevImg, InputArray nextImg,InputOutputArray flow, double pyrScale, int levels, int winsize, intiterations, int polyN, double polySigma, int flags)
大部分参数在论文中都有一套比较好的值的，直接采用他们的就好了
// Farneback dense optical flow calculate and show in Munsell system of colors
// Author : Zouxy
// Date : 2013-3-15
// HomePage : http://blog.csdn.net/zouxy09
// Email : zouxy09@qq.com
// API calcOpticalFlowFarneback() comes from OpenCV, and this
// 2D dense optical flow algorithm from the following paper:
// Gunnar Farneback. "Two-Frame Motion Estimation Based on Polynomial Expansion".
// And the OpenCV source code locate in ..opencv2.4.3modulesvideosrcoptflowgf.cpp
#include <iostream>
#include "opencv2/opencv.hpp"
using namespace cv;
using namespace std;
#define UNKNOWN_FLOW_THRESH 1e9
// Color encoding of flow vectors from:
// http://members.shaw.ca/quadibloc/other/colint.htm
// This code is modified from:
// http://vision.middlebury.edu/flow/data/
void makecolorwheel(vector<Scalar> &colorwheel)
{
int RY = 15;
int YG = 6;
int GC = 4;
int CB = 11;
int BM = 13;
int MR = 6;
int i;
for (i = 0; i < RY; i++) colorwheel.push_back(Scalar(255, 255*i/RY, 0));
for (i = 0; i < YG; i++) colorwheel.push_back(Scalar(255-255*i/YG, 255, 0));
for (i = 0; i < GC; i++) colorwheel.push_back(Scalar(0, 255, 255*i/GC));
for (i = 0; i < CB; i++) colorwheel.push_back(Scalar(0, 255-255*i/CB, 255));
for (i = 0; i < BM; i++) colorwheel.push_back(Scalar(255*i/BM, 0, 255));
for (i = 0; i < MR; i++) colorwheel.push_back(Scalar(255, 0, 255-255*i/MR));
}
void motionToColor(Mat flow, Mat &color)
{
if (color.empty())
color.create(flow.rows, flow.cols, CV_8UC3);
static vector<Scalar> colorwheel; //Scalar r,g,b
if (colorwheel.empty())
makecolorwheel(colorwheel);
// determine motion range:
float maxrad = -1;
// Find max flow to normalize fx and fy
for (int i= 0; i < flow.rows; ++i)
{
for (int j = 0; j < flow.cols; ++j)
{
Vec2f flow_at_point = flow.at<Vec2f>(i, j);
float fx = flow_at_point[0];
float fy = flow_at_point[1];
if ((fabs(fx) > UNKNOWN_FLOW_THRESH) || (fabs(fy) > UNKNOWN_FLOW_THRESH))
continue;
float rad = sqrt(fx * fx + fy * fy);
maxrad = maxrad > rad ? maxrad : rad;
}
}
for (int i= 0; i < flow.rows; ++i)
{
for (int j = 0; j < flow.cols; ++j)
{
uchar *data = color.data + color.step[0] * i + color.step[1] * j;
Vec2f flow_at_point = flow.at<Vec2f>(i, j);
float fx = flow_at_point[0] / maxrad;
float fy = flow_at_point[1] / maxrad;
if ((fabs(fx) > UNKNOWN_FLOW_THRESH) || (fabs(fy) > UNKNOWN_FLOW_THRESH))
{
data[0] = data[1] = data[2] = 0;
continue;
}
float rad = sqrt(fx * fx + fy * fy);
float angle = atan2(-fy, -fx) / CV_PI;
float fk = (angle + 1.0) / 2.0 * (colorwheel.size()-1);
int k0 = (int)fk;
int k1 = (k0 + 1) % colorwheel.size();
float f = fk - k0;
//f = 0; // uncomment to see original color wheel
for (int b = 0; b < 3; b++)
{
float col0 = colorwheel[k0][b] / 255.0;
float col1 = colorwheel[k1][b] / 255.0;
float col = (1 - f) * col0 + f * col1;
if (rad <= 1)
col = 1 - rad * (1 - col); // increase saturation with radius
else
col *= .75; // out of range
data[2 - b] = (int)(255.0 * col);
}
}
}
}
int main(int, char**)
{
VideoCapture cap;
cap.open(0);
//cap.open("test_02.wmv");
if( !cap.isOpened() )
return -1;
Mat prevgray, gray, flow, cflow, frame;
namedWindow("flow", 1);
Mat motion2color;
for(;;)
{
double t = (double)cvGetTickCount();
cap >> frame;
cvtColor(frame, gray, CV_BGR2GRAY);
imshow("original", frame);
if( prevgray.data )
{
calcOpticalFlowFarneback(prevgray, gray, flow, 0.5, 3, 15, 3, 5, 1.2, 0);
motionToColor(flow, motion2color);
imshow("flow", motion2color);
}
if(waitKey(10)>=0)
break;
std::swap(prevgray, gray);
t = (double)cvGetTickCount() - t;
cout << "cost time: " << t / ((double)cvGetTickFrequency()*1000.) << endl;
}
return 0;
}
// Farneback dense optical flow calculate and show in Munsell system of colors
// Author : Zouxy
// Date : 2013-3-15
// HomePage : http://blog.csdn.net/zouxy09
// Email : zouxy09@qq.com
// API calcOpticalFlowFarneback() comes from OpenCV, and this
// 2D dense optical flow algorithm from the following paper:
// Gunnar Farneback. "Two-Frame Motion Estimation Based on Polynomial Expansion".
// And the OpenCV source code locate in ..opencv2.4.3modulesvideosrcoptflowgf.cpp
#include <iostream>
#include "opencv2/opencv.hpp"
using namespace cv;
using namespace std;
#define UNKNOWN_FLOW_THRESH 1e9
// Color encoding of flow vectors from:
// http://members.shaw.ca/quadibloc/other/colint.htm
// This code is modified from:
// http://vision.middlebury.edu/flow/data/
void makecolorwheel(vector<Scalar> &colorwheel)
{
int RY = 15;
int YG = 6;
int GC = 4;
int CB = 11;
int BM = 13;
int MR = 6;
int i;
for (i = 0; i < RY; i++) colorwheel.push_back(Scalar(255, 255*i/RY, 0));
for (i = 0; i < YG; i++) colorwheel.push_back(Scalar(255-255*i/YG, 255, 0));
for (i = 0; i < GC; i++) colorwheel.push_back(Scalar(0, 255, 255*i/GC));
for (i = 0; i < CB; i++) colorwheel.push_back(Scalar(0, 255-255*i/CB, 255));
for (i = 0; i < BM; i++) colorwheel.push_back(Scalar(255*i/BM, 0, 255));
for (i = 0; i < MR; i++) colorwheel.push_back(Scalar(255, 0, 255-255*i/MR));
}
void motionToColor(Mat flow, Mat &color)
{
if (color.empty())
color.create(flow.rows, flow.cols, CV_8UC3);
static vector<Scalar> colorwheel; //Scalar r,g,b
if (colorwheel.empty())
makecolorwheel(colorwheel);
// determine motion range:
float maxrad = -1;
// Find max flow to normalize fx and fy
for (int i= 0; i < flow.rows; ++i)
{
for (int j = 0; j < flow.cols; ++j)
{
Vec2f flow_at_point = flow.at<Vec2f>(i, j);
float fx = flow_at_point[0];
float fy = flow_at_point[1];
if ((fabs(fx) > UNKNOWN_FLOW_THRESH) || (fabs(fy) > UNKNOWN_FLOW_THRESH))
continue;
float rad = sqrt(fx * fx + fy * fy);
maxrad = maxrad > rad ? maxrad : rad;
}
}
for (int i= 0; i < flow.rows; ++i)
{
for (int j = 0; j < flow.cols; ++j)
{
uchar *data = color.data + color.step[0] * i + color.step[1] * j;
Vec2f flow_at_point = flow.at<Vec2f>(i, j);
float fx = flow_at_point[0] / maxrad;
float fy = flow_at_point[1] / maxrad;
if ((fabs(fx) > UNKNOWN_FLOW_THRESH) || (fabs(fy) > UNKNOWN_FLOW_THRESH))
{
data[0] = data[1] = data[2] = 0;
continue;
}
float rad = sqrt(fx * fx + fy * fy);
float angle = atan2(-fy, -fx) / CV_PI;
float fk = (angle + 1.0) / 2.0 * (colorwheel.size()-1);
int k0 = (int)fk;
int k1 = (k0 + 1) % colorwheel.size();
float f = fk - k0;
//f = 0; // uncomment to see original color wheel
for (int b = 0; b < 3; b++)
{
float col0 = colorwheel[k0][b] / 255.0;
float col1 = colorwheel[k1][b] / 255.0;
float col = (1 - f) * col0 + f * col1;
if (rad <= 1)
col = 1 - rad * (1 - col); // increase saturation with radius
else
col *= .75; // out of range
data[2 - b] = (int)(255.0 * col);
}
}
}
}
int main(int, char**)
{
VideoCapture cap;
cap.open(0);
//cap.open("test_02.wmv");
if( !cap.isOpened() )
return -1;
Mat prevgray, gray, flow, cflow, frame;
namedWindow("flow", 1);
Mat motion2color;
for(;;)
{
double t = (double)cvGetTickCount();
cap >> frame;
cvtColor(frame, gray, CV_BGR2GRAY);
imshow("original", frame);
if( prevgray.data )
{
calcOpticalFlowFarneback(prevgray, gray, flow, 0.5, 3, 15, 3, 5, 1.2, 0);
motionToColor(flow, motion2color);
imshow("flow", motion2color);
}
if(waitKey(10)>=0)
break;
std::swap(prevgray, gray);
t = (double)cvGetTickCount() - t;
cout << "cost time: " << t / ((double)cvGetTickFrequency()*1000.) << endl;
}
return 0;
}
这个是效果：
一个挥动的手：
虽然也有背景在动，但是因为他们的运动方向不一样，所以还是可以辨认出来前面那个是手，在前景和背景运动不统一的时候，还是可以辨认出来的。
