光流算法从理论到实践专题1_cv::cuda::sparsepyrlkopticalflow-优快云博客

资料搜索

本人总结

1、什么是光流？

通俗来说就是，在物体运动过程中，会造成对应位置的亮度变化。例如，在停车场中地面为灰色的某一位置，停放着一辆红色保时捷，我们看该位置为红色。但过了一分钟后这辆保时捷开走了，该位置颜色变为灰色，我们在看是该位置为灰色。说白了就是，反映在我们人眼的是颜色、亮度发生变化，这就是光流。从图像角度来讲，视频图像的一帧中代表同一对象的像素点移到下一帧时该对象的像素移动量。

2、为什么要研究光流？

光流概念是Gibson在1950年首先提出来的，发展到很多研究院对该算法进行了各种的改进，其中最新的光流算法，就是近几年通过引入神经网络提出的FlowNet/FlowNet2。为什么会到现在还是有很多科研工作者对光流算法进行研究呢！我们可以举个例子来说明这个问题：假如我们在漆黑的夜晚行走时，我们看见远方来了一个正向你逐渐靠近的人，因为周围漆黑你无法判断这个人是谁，但是当你观察几秒，发现走路的姿态可以初步判定这个人是赵三。当这个人近距离走进你的时候，发现这个人确实是赵三。从上面这个例子可以分析到，我们只考虑单张图像去做分析，有的时候很难做出结论，但是当我们考虑时间相关的一段视频，更有助于我们做分析和判断。而且随着5G技术的发展，视频传输不在是问题，计算机视觉从单张图像研究也会逐渐向视频研究展开，光流估计作为视频理解的就显得重中之重了。

3、光流算法推导及实现

3.1 Lucas-Kanade

LK光流算法是在1981年提出，想看它对应论文的可以点击标题，就可以跳转到对应的pdf。LK算法是基于三个假设（亮度恒定、运动是小运动、空间一致）展开研究的，我们下面主要分析这三个假设的数学推导，以及思考一下为什么会使这三点。

3.1.1 亮度恒定

（1）概念

说白了就是我们获得的视频帧率非常快，按照一般情况下视频帧率为25帧，所以两帧之间的时间间隔也就40ms，而我们可以大体认为40ms同一位置亮度的变化是基本不变的。

（2）公式转化

根据上述的描述，我们可以将其转化为如下数学公式： $I(x,y,t)=I(x+u,y+v,t+\Delta t)$ ，即t时刻（x,y）位置的亮度与 $t+\Delta t$ 时刻对应位置 $(x+u,y+v)$ 的亮度相同。

上面的假设我们可以得到一个等式方程，但是我们要如何求解这个方程呢。因为一般视频中两帧之间的间隔非常短，即认为两帧之间运动非常小。

3.1.2 运动是小运动

（1）概念

即运动的变化可以理解为是时间的导数

（2）公式展开

根据泰勒展开可以得到 $I(x+u,y+v,t+\Delta t) = I(x,y,t) + I_{x}^{'}u + I_{y}^{'}v + I_{t}^{'} \Delta t$ 。关于泰勒展开的相关知识，可以看此链接。

根据3.1.1中的等式可以将公式转化为： $I(x,y,t) = I(x,y,t) + I_{x}^{'}u + I_{y}^{'}v + I_{t}^{'} \Delta t$

由此可以得到公式： $I_{x}^{'}u + I_{y}^{'}v + I_{t}^{'} \Delta t = 0$

转化为矩阵表达： $\begin{bmatrix} I_{x}^{'} & I_{y}^{'} \end{bmatrix}\begin{bmatrix} v\\ u \end{bmatrix}= -I_{t}^{'} \Delta t$

其中， $I_{x}^{'}$ 和 $I_{y}^{'}$ 分别对应在x方向和y方向的偏导数， $I_{t}^{'}$ 对应在t时刻在（x，y）处像素的时间导数。 $I_{t}^{'}\Delta t$

在位置（x,y）的亮度差，即可以表示为 $\Delta I_{t}^{}$ 。

因此公式可以转化为： $\begin{bmatrix} I_{x}^{'} & I_{y}^{'} \end{bmatrix}\begin{bmatrix} v\\ u \end{bmatrix}= -\Delta I_{t}^{}$

给定两张图像，我们可以已知的是 $I_{x}^{'}$ 、 $I_{y}^{'}$ 、 $\Delta I_{t}^{}$ 。带求解的是u、v值，但是一个位置点只能得到一个方程，所以如何求解位置数u、v。我们肯定都会想到就是增加等式方程，但我们要怎么增加等式方程，我们都知道物以类聚，所以只要距离比较小的，我们就会认为这类性质一般都是一样的。所以，我们又提出下面的假设：空间一致

3.1.3 空间一致

（1）概念

根据上面3.1.2中叙述的一样，取距离最近的点认为它们的性质是一致的，即满足同一个方程。因为我们可以认为（x,y）周围3*3大小或者其它大小空间中的周围点满足同一个方程。

（2）公式

根据概念可以转化为下面的公式，其中1,2,...,n代表是（x,y）周围的像素点。

$\begin{bmatrix} I_{x}^{'(1)}& I_{y}^{'(1)}\\ I_{x}^{'(2)}& I_{y}^{'(2)} \\ .....&..... \\ I_{x}^{'(n)}& I_{y}^{'(n)} \end{bmatrix} \begin{bmatrix} v\\ u \end{bmatrix}= \begin{bmatrix} -\Delta I_{t}^{(1)}\\ -\Delta I_{t}^{(1)}\\ .....\\ -\Delta I_{t}^{(n)}\\ \end{bmatrix}$

要求解这个超参数方程，一般做法是使用最小二乘法求解，不了解最小二乘法的看此链接。因此，上面方程转化为下面形式：

$Ax=b$ ,其中A代表 $\begin{bmatrix} I_{x}^{'(1)}& I_{y}^{'(1)}\\ I_{x}^{'(2)}& I_{y}^{'(2)} \\ .....&..... \\ I_{x}^{'(n)}& I_{y}^{'(n)} \end{bmatrix}$ ，x代表 $\begin{bmatrix} v\\ u \end{bmatrix}$ ，b代表 $\begin{bmatrix} -\Delta I_{t}^{(1)}\\ -\Delta I_{t}^{(1)}\\ .....\\ -\Delta I_{t}^{(n)}\\ \end{bmatrix}$

根据线性代数矩阵运算知识（不懂得看此链接）可得 $x=(A^{T}A)^{-1}A^{T}b$

为了使方程有解， $(A^{T}A)^{}$ 要求是可逆的。但是当 $(A^{T}A)^{}$ 不可逆时，即可能会出现孔径的问题，如下图所示情况。所以在 $(A^{T}A)^{}$ 可逆的情况下，一般都是亮度变化比较明显的位置。我们通常所用的Harris角点检测（不懂得看此链接，后面有时间我也会单独写一篇角点检测专题），就是用的是 $(A^{T}A)^{}$ 可逆的性质。

3.2 OpenCV实现

下面是使用OpenCV例子，该程序中包含稠密光流和稀疏光流，只需要flow_type即可，测试了针对当前这个场景效果差不多。

#include <iostream>
#include <vector>

#include <opencv2/core.hpp>
#include <opencv2/core/utility.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/video.hpp>
#include <opencv2/cudaoptflow.hpp>
#include <opencv2/cudaimgproc.hpp>
#include <opencv2/cudaarithm.hpp>

using namespace std;
using namespace cv;
using namespace cv::cuda;

static void download(const GpuMat& d_mat, vector<Point2f>& vec)
{
	vec.resize(d_mat.cols);
	Mat mat(1, d_mat.cols, CV_32FC2, (void*)&vec[0]);
	d_mat.download(mat);
}

static void download(const GpuMat& d_mat, vector<uchar>& vec)
{
	vec.resize(d_mat.cols);
	Mat mat(1, d_mat.cols, CV_8UC1, (void*)&vec[0]);
	d_mat.download(mat);
}

static void drawArrows(Mat& frame, const vector<Point2f>& prevPts, const vector<Point2f>& nextPts, const vector<uchar>& status, Scalar line_color = Scalar(0, 0, 255))
{
	for (size_t i = 0; i < prevPts.size(); ++i)
	{
		if (status[i])
		{
			int line_thickness = 1;

			Point p = prevPts[i];
			Point q = nextPts[i];

			double angle = atan2((double)p.y - q.y, (double)p.x - q.x);

			double hypotenuse = sqrt((double)(p.y - q.y)*(p.y - q.y) + (double)(p.x - q.x)*(p.x - q.x));

			if (hypotenuse < 1.0)
				continue;

			// Here we lengthen the arrow by a factor of three.
			q.x = (int)(p.x - 3 * hypotenuse * cos(angle));
			q.y = (int)(p.y - 3 * hypotenuse * sin(angle));

			// Now we draw the main line of the arrow.
			line(frame, p, q, line_color, line_thickness);

			// Now draw the tips of the arrow. I do some scaling so that the
			// tips look proportional to the main line of the arrow.

			p.x = (int)(q.x + 9 * cos(angle + CV_PI / 4));
			p.y = (int)(q.y + 9 * sin(angle + CV_PI / 4));
			line(frame, p, q, line_color, line_thickness);

			p.x = (int)(q.x + 9 * cos(angle - CV_PI / 4));
			p.y = (int)(q.y + 9 * sin(angle - CV_PI / 4));
			line(frame, p, q, line_color, line_thickness);
		}
	}
}

inline bool isFlowCorrect(Point2f u)
{
	return !cvIsNaN(u.x) && !cvIsNaN(u.y) && fabs(u.x) < 1e9 && fabs(u.y) < 1e9;
}

static Vec3b computeColor(float fx, float fy)
{
	static bool first = true;

	// relative lengths of color transitions:
	// these are chosen based on perceptual similarity
	// (e.g. one can distinguish more shades between red and yellow
	//  than between yellow and green)
	const int RY = 15;
	const int YG = 6;
	const int GC = 4;
	const int CB = 11;
	const int BM = 13;
	const int MR = 6;
	const int NCOLS = RY + YG + GC + CB + BM + MR;
	static Vec3i colorWheel[NCOLS];

	if (first)
	{
		int k = 0;

		for (int i = 0; i < RY; ++i, ++k)
			colorWheel[k] = Vec3i(255, 255 * i / RY, 0);

		for (int i = 0; i < YG; ++i, ++k)
			colorWheel[k] = Vec3i(255 - 255 * i / YG, 255, 0);

		for (int i = 0; i < GC; ++i, ++k)
			colorWheel[k] = Vec3i(0, 255, 255 * i / GC);

		for (int i = 0; i < CB; ++i, ++k)
			colorWheel[k] = Vec3i(0, 255 - 255 * i / CB, 255);

		for (int i = 0; i < BM; ++i, ++k)
			colorWheel[k] = Vec3i(255 * i / BM, 0, 255);

		for (int i = 0; i < MR; ++i, ++k)
			colorWheel[k] = Vec3i(255, 0, 255 - 255 * i / MR);

		first = false;
	}

	const float rad = sqrt(fx * fx + fy * fy);
	const float a = atan2(-fy, -fx) / (float)CV_PI;

	const float fk = (a + 1.0f) / 2.0f * (NCOLS - 1);
	const int k0 = static_cast<int>(fk);
	const int k1 = (k0 + 1) % NCOLS;
	const float f = fk - k0;

	Vec3b pix;

	for (int b = 0; b < 3; b++)
	{
		const float col0 = colorWheel[k0][b] / 255.0f;
		const float col1 = colorWheel[k1][b] / 255.0f;

		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

		pix[2 - b] = static_cast<uchar>(255.0 * col);
	}

	return pix;
}

static void drawOpticalFlow(const Mat_<float>& flowx, const Mat_<float>& flowy, Mat& dst, float maxmotion = -1)
{
	dst.create(flowx.size(), CV_8UC3);
	dst.setTo(Scalar::all(0));

	// determine motion range:
	float maxrad = maxmotion;

	if (maxmotion <= 0)
	{
		maxrad = 1;
		for (int y = 0; y < flowx.rows; ++y)
		{
			for (int x = 0; x < flowx.cols; ++x)
			{
				Point2f u(flowx(y, x), flowy(y, x));

				if (!isFlowCorrect(u))
					continue;

				maxrad = max(maxrad, sqrt(u.x * u.x + u.y * u.y));
			}
		}
	}

	for (int y = 0; y < flowx.rows; ++y)
	{
		for (int x = 0; x < flowx.cols; ++x)
		{
			Point2f u(flowx(y, x), flowy(y, x));

			if (isFlowCorrect(u))
				dst.at<Vec3b>(y, x) = computeColor(u.x / maxrad, u.y / maxrad);
		}
	}
}

static void showFlow(const char* name, const GpuMat& d_flow)
{
	GpuMat planes[2];
	cuda::split(d_flow, planes);

	Mat flowx(planes[0]);
	Mat flowy(planes[1]);

	Mat out;
	drawOpticalFlow(flowx, flowy, out, 10);

	imshow(name, out);
}

template <typename T> inline T clamp(T x, T a, T b)
{
	return ((x) > (a) ? ((x) < (b) ? (x) : (b)) : (a));
}

template <typename T> inline T mapValue(T x, T a, T b, T c, T d)
{
	x = clamp(x, a, b);
	return c + (d - c) * (x - a) / (b - a);
}

int main(int argc, const char* argv[])
{
	const char* keys =
		"{ h             help   |        | print help message }"
		"{ l             left   | ../data/pic1.png       | specify left image }"
		"{ r             right  | ../data/pic2.png       | specify right image }"
		"{ flow                 | sparse | specify flow type [PyrLK] }"
		"{ gray                 |        | use grayscale sources [PyrLK Sparse] }"
		"{ win_size             | 21     | specify windows size [PyrLK] }"
		"{ max_level            | 3      | specify max level [PyrLK] }"
		"{ iters                | 30     | specify iterations count [PyrLK] }"
		"{ points               | 4000   | specify points count [GoodFeatureToTrack] }"
		"{ min_dist             | 0      | specify minimal distance between points [GoodFeatureToTrack] }";

	CommandLineParser cmd(argc, argv, keys);

	if (cmd.has("help") || !cmd.check())
	{
		cmd.printMessage();
		cmd.printErrors();
		return 0;
	}

	string fname0 = "./basketball1.png";//cmd.get<string>("left");
	string fname1 = "./basketball2.png";//cmd.get<string>("right");


	string flow_type = "dense";//cmd.get<string>("flow");
	bool is_sparse = true;
	if (flow_type == "sparse")
	{
		is_sparse = true;
	}
	else if (flow_type == "dense")
	{
		is_sparse = false;
	}
	else
	{
		cerr << "please specify 'sparse' or 'dense' as flow type" << endl;
		return -1;
	}

	bool useGray = 0;//cmd.has("gray");
	int winSize = 21;//cmd.get<int>("win_size");
	int maxLevel = 3;//cmd.get<int>("max_level");
	int iters = 30;//cmd.get<int>("iters");
	int points = 4000;//cmd.get<int>("points");
	double minDist = 0;//cmd.get<double>("min_dist");

	Mat frame0 = imread(fname0);
	Mat frame1 = imread(fname1);

	if (frame0.empty() || frame1.empty())
	{
		cout << "Can't load input images" << endl;
		return -1;
	}

	cout << "Image size : " << frame0.cols << " x " << frame0.rows << endl;
	cout << "Points count : " << points << endl;

	cout << endl;

	Mat frame0Gray;
	cv::cvtColor(frame0, frame0Gray, COLOR_BGR2GRAY);
	Mat frame1Gray;
	cv::cvtColor(frame1, frame1Gray, COLOR_BGR2GRAY);

	// goodFeaturesToTrack
	GpuMat d_frame0Gray(frame0Gray);
	GpuMat d_prevPts;

	Ptr<cuda::CornersDetector> detector = cuda::createGoodFeaturesToTrackDetector(d_frame0Gray.type(), points, 0.01, minDist);
	detector->detect(d_frame0Gray, d_prevPts);

	GpuMat d_frame0(frame0);
	GpuMat d_frame1(frame1);
	GpuMat d_frame1Gray(frame1Gray);
	GpuMat d_nextPts;
	GpuMat d_status;
	GpuMat d_flow(frame0.size(), CV_32FC2);

	if (is_sparse)
	{
		// Sparse
		Ptr<cuda::SparsePyrLKOpticalFlow> d_pyrLK_sparse = cuda::SparsePyrLKOpticalFlow::create(
			Size(winSize, winSize), maxLevel, iters);
		d_pyrLK_sparse->calc(useGray ? d_frame0Gray : d_frame0, useGray ? d_frame1Gray : d_frame1, d_prevPts, d_nextPts, d_status);

		// Draw arrows
		vector<Point2f> prevPts(d_prevPts.cols);
		download(d_prevPts, prevPts);

		vector<Point2f> nextPts(d_nextPts.cols);
		download(d_nextPts, nextPts);

		vector<uchar> status(d_status.cols);
		download(d_status, status);

		namedWindow("PyrLK [Sparse]", WINDOW_NORMAL);
		drawArrows(frame0, prevPts, nextPts, status, Scalar(255, 0, 0));
		imshow("PyrLK [Sparse]", frame0);
	}
	else
	{
		// Dense
		Ptr<cuda::DensePyrLKOpticalFlow> d_pyrLK_dense = cuda::DensePyrLKOpticalFlow::create(
			Size(winSize, winSize), maxLevel, iters);
		d_pyrLK_dense->calc(d_frame0Gray, d_frame1Gray, d_flow);

		// Draw flows
		namedWindow("PyrLK [Dense] Flow Field", WINDOW_NORMAL);
		showFlow("PyrLK [Dense] Flow Field", d_flow);
	}

	waitKey(0);

	return 0;
}

3.3 实现效果

4、常见的光流

（1）基于梯度方法

基于梯度的方法又称为微分法，利用时变图像灰度（或者其滤波后的形式）进行时空微分（时空梯度函数）来计算像素的速度矢量。其中，最具代表的就是Horn-Schunck算法与Lucas-Kanade(LK)算法。Horn-Schunck算法在光流约束方程的基础上附加了全局平滑假设，假设整张图上光流是光滑，即物体运动矢量是平滑或者是缓慢变化的。后面很多人在此基础上进行改进，Nage采用有条件的平滑约束，即通过加权矩阵的控制对梯度进行不同平滑处理；Black和Anandan针对多运动提出分段平滑方法。

优点：计算简单和可以达到不错的效果

（2）基于匹配的方法

该方法目前主要有两种方法：基于特征和区域。其中基于特征方法主要是对目标的主要特征进行定位和跟踪，对大目标和特征明显具有很强的鲁棒性；基于区域方法主要是对类似的区域进行定位，然后通过相似区域的位移进行光流处理，主要在视频编码中使用；两者方法都是因为稀疏所以特征提取和匹配都很难，即使是亚像素的计算量就很大。

（3）基于能量的方法

该方法又称为基于频率的方法，对图像进行时空滤波处理，即对时间和空间进行整合；

缺点：降低了光流时间和空间分辨率；可靠性评价比较难；计算复杂度很大；

（4）基于相位的方法

Fleet和Jepson使用相位的信息来进行光流计算处理；

优点：对图像序列适用范围很宽，而且速度估计比较精确；

缺点：模型设计合理，但是时间复杂度很高；对图像序列的时间混叠比较敏感；

（5）神经动力学方法

该方法使用神经网络建立视觉运动感知的神经动力模型，让它对生物视觉系统功能与结构比较直接的模拟。

更多《计算机视觉与图形学》知识，可关注下方公众号：