C++实现LK光流法/计算稀疏光流

我又来了,这万恶的大作业!

1 程序简介

自行实现LK光流计算,实现了对指定视频(三通道彩色)计算稀疏光流并显示在图像上。
主要包括的自编函数有:
RGB2GRAY (),将彩色图转为灰度图。
ScharrConvolution (),使用Scharr算子求XY方向梯度。
getR(),获得角点响应值R,用于获取最优角点。
NMS(),非极大值抑制。
myCompare(),自定义排序函数。
getBestCorner(),根据R值由大到小对坐标进行排序,获取前n个R值最大的坐标,若整体不足n个则全部返回。
getGradT(),计算t方向梯度。
getUV(),计算UV速度。

2 算法流程

首先计算图像在X、Y和T方向的梯度,然后计算角点响应值R,进行最大值抑制,选出最优的n个角点,到这里就计算出了全部的角点,之后不再更新。由角点坐标和三向梯度根据公式计算各点的速度。
进入循环,先取出下一帧图片,由坐标和速度在其上画出光流,然后重复上方计算速度的步骤算出新帧的速度,之后再取出下一帧图片,重复这些操作。

2.1 彩色图转灰度图

以0.114:0.587:0.299的系数对BGR加权求和。

2.2 计算X、Y方向梯度

分别使用ScharrX和ScharrY算子对图像进行卷积即可获得两方向的梯度。

2.3 计算角点响应值R

根据X向和Y向梯度计算出各坐标的角点响应值。

2.4 非极大值抑制

简称NMS,会剔除一定大小的区域内的非极大值。

2.5 获取最优角点

函数名为getBestCorner(),传入角点序列和期望获取的最优角点数。函数会根据坐标对应的R值从大到小对坐标点进行排序,取前n个角点,若整体不足或等于n个则全部返回。

void getBestCorner(vector<vector<double>>& corner_list, int n) {
	if (corner_list.size() <= n) {  //若角点数量不足或等于n,则退出函数
		return;
	}
	sort(corner_list.begin(), corner_list.end(), myCompare);
	corner_list.erase(corner_list.begin() + n, corner_list.end());
	//for (int i = 0; i < corner_list.size(); i++) {  //用于查看所有坐标,可在选中后,先后按ctrl+K和ctrl+U来取消注释(Visual Studio)
	//	cout << "x:" << corner_list[i][0] << ",y:" << corner_list[i][1] << ",R:" << corner_list[i][2] << endl;
	//}
	cout << "已提取前" << corner_list.size() << "个角点,目标个数为:" << n << endl;

	return;
}

2.6 计算t方向梯度

实现较为简单,只需计算前后两帧图像坐标的差的绝对值即可。

void getGradT(Mat img, Mat img_next, float* grad_T) {
	uchar* p_img = img.data;
	uchar* p_img_next = img_next.data;
	for (int i = 0; i < frame_height; i++) {//卷积
		for (int j = 0; j < frame_width; j++) {
			grad_T[i * frame_width + j] = abs(p_img_next[i * frame_width + j] - p_img[i * frame_width + j]);
		}
	}
}

2.7 计算坐标点速度

算法中最重要的函数,根据传入的三向梯度(X、Y和T)计算出坐标点对应的速度。
函数首先剔除边缘点,即即将运动出图像范围的角点,防止指针指向非法地址。

if (row + window / 2 > frame_height - 1 || row - window / 2 < 0 || col + window / 2 > frame_width - 1 || col - window / 2 < 0) {
	corner_list.erase(corner_list.begin() + i);  //如果清除了第i个元素,则后面的元素会递补第i个位置,所以无需i++
}
else {
	i++;
}

然后根据如下公式(最小二乘法等式)计算速度:

如今已知Ix、Iy和It,只需计算出左侧22矩阵的逆矩阵即可计算出速度值,而22矩阵的逆矩阵易得:

故可以计算出速度u、v。

for (int m = -window / 2; m < window / 2 + 1; m++) {
	for (int n = -window / 2; n < window / 2 + 1; n++) {
		gradXX += grad_X[(row + m) * frame_width + col + n] * grad_X[(row + m) * frame_width + col + n];
		gradXY += grad_X[(row + m) * frame_width + col + n] * grad_Y[(row + m) * frame_width + col + n];
		gradYY += grad_Y[(row + m) * frame_width + col + n] * grad_Y[(row + m) * frame_width + col + n];
		gradXT += grad_X[(row + m) * frame_width + col + n] * grad_T[(row + m) * frame_width + col + n];
		gradYT += grad_Y[(row + m) * frame_width + col + n] * grad_T[(row + m) * frame_width + col + n];
	}
}
double deno = gradXX * gradYY - gradXY * gradXY;  //求矩阵逆生成的分母
double U = (gradYY * gradXT + (-gradXY) * gradYT) / deno;
double V = ((-gradXY) * gradXT + gradXX * gradYT) / deno;

3 实现细节

3.1 获得最优角点细节

当图像过于平滑时,角点数可能会很少,如果还选取前n个就可能会因为超出变量范围导致报错,所以需要在函数开始时加入if语句判断角点数量。

if (corner_list.size() <= n) {  //若角点数量不足或等于n,则退出函数
		return;
}

3.2 计算坐标点速度细节

计算速度需要卷积,而卷积最可能发生异常的地方就是索引超出范围,即点跑到图像外,那么我们在卷积前就需要对每个角点逐一判断在卷积时是否会超出边界,这个判断过程实现较为简单,但有一个细节,即索引i的变化。假如我们判断到第6个角点,即i=5,我们发现这个角点在卷积时会出界,所以我们调用erase函数删除这个值,而在删除这个值之后,原处于第7个位置的坐标补到第6个位置,如果我们此时按照常规仍然对i进行加一操作,那么此时i=6,我们就会略过此时处于i=5位置处的坐标,没有判断它在卷积时是否会出界,那么我们的程序在卷积时就有可能发生内存访问异常的报错。
要解决这个问题也很简单,即将i++改为仅当判断坐标不会出界后执行,这样当坐标被删除时就不会执行i++,我们也就不会略过任何一个坐标点。

for (int i = 0;;) {
	row = corner_list[i][0] + 0.5;
	col = corner_list[i][1] + 0.5;
	if (row + window / 2 > frame_height - 1 || row - window / 2 < 0 || col + window / 2 > frame_width - 1 || col - window / 2 < 0) {
		corner_list.erase(corner_list.begin() + i);  //如果清除了第i个元素,则后面的元素会递补第i个位置,所以无需i++
	}
	else {
		i++;
	}
	if (i >= corner_list.size()) {
		break;
	}
}

后面计算速度只需要对照公式实现计算过程,并没有值得一提的细节。

3.3 主函数细节

首先是由于视频尺寸过大会导致处理速度过慢、变量超出允许范围等问题,所以我们在读入每张图片时都需要减小尺寸。
Size frame_size = Size(870, 540);
其次是如何储存点坐标。这里可以使用vector容器,也可以使用float指针,不过如果想使用float数组指针就必须创建动态数组(因为数组长度未定,是根据卷积核大小来确定的),还要注意释放的问题。这里选择了使用vector容器,编程比较方便,可以自由增删元素,这样就不必为了防止内存访问异常而定义一个很大的数组。
vector<vector> corner_list;//储存角点坐标
然后是NMS窗口大小的确定。原则上来说是根据图像内物体密度和图像大小来综合确定的,这里经过测试选择了11,效果比较好。
NMS(result_img, corner_list, 11);
还有比较难的一点是如何在图像上画出光流。我这里采用的方法为:定义一个名为mask的Mat变量,初始化全黑,即像素全为0,然后就在mask上绘制光流,每次更新光流图像时只需要将mask与实时图像相加即可。但由于进行的是相加操作,我们就不能自定义光流颜色,只能将光流设为全白,即像素值全为255,这样在mask与视频图像相加时,光流对应的区域的像素值也会被设置为255(因为像素值上限就为255)。
另外在测试过程中发现光流很容易由于速度波动导致跟不上实际物体,所以我在由速度更新角点坐标时将速度乘以二。因为坐标点超过了实际位置不会导致跟踪目标丢失,他会留在原地等待实际物体运动到坐标点处,但如果坐标点落后于物体,那这个点只能停留在原地了。

corner_list[i][0] += 2 * velocity_list[i][0];
corner_list[i][1] += 2 * velocity_list[i][1];

最后释放动态数组。

4 源代码

#include <opencv2/opencv.hpp>
#include <opencv2/core/utils/logger.hpp>  //用来关闭opencv的日志输出的(日志看着比较乱)

using namespace cv;
using namespace std;

int frame_count, frame_fps, frame_duration, frame_height, frame_width, frame_channels = 3;  //储存图像帧数、帧率、高度、宽度、通道数

void RGB2GRAY(Mat& img, Mat& GRAY_img);

void ScharrConvolution(Mat& img, float* grad_X, float* grad_Y);

void getR(double* result_img, float* grad_X, float* grad_Y);

void NMS(double* result_img, vector<vector<double>>& corner_list, int core = 40);

bool myCompare(vector<double>& corner_list1, vector<double>& corner_list2);

void getBestCorner(vector<vector<double>>& corner_list, int n);

void getGradT(Mat img, Mat img_next, float* grad_T);

void getUV(float* grad_X, float* grad_Y, float* grad_T, vector<vector<double>>& corner_list, vector<vector<double>>& velcity_list, int window = 7);


int main(){
	cv::utils::logging::setLogLevel(utils::logging::LOG_LEVEL_SILENT);  //关闭日志输出

	//获取视频信息
	VideoCapture capture("C:\\Users\\33507\\Desktop\\123.mp4");
	frame_count = int(capture.get(CAP_PROP_FRAME_COUNT));
	frame_fps = int(capture.get(CAP_PROP_FPS));
	frame_height = int(capture.get(CAP_PROP_FRAME_HEIGHT));
	frame_width = int(capture.get(CAP_PROP_FRAME_WIDTH));
	cout << "时长:" << frame_count / frame_fps << "s,帧率:" << frame_fps << "Hz,原宽高:" << frame_width << "," << frame_height << endl;

	//由于视频尺寸过大会导致处理速度过慢、变量超出允许范围等问题,故减小视频尺寸
	frame_height = 540;
	frame_width = 870;
	Size frame_size = Size(870, 540);

	//变量初始化
	vector<vector<double>> corner_list;//储存角点坐标
	vector<vector<double>> velocity_list;//储存角点速度
	Mat frame, frame_next, img, img_next, img_GRAY, img_GRAY_next, mask;
	double* result_img = new double[frame_width * frame_height];
	float* grad_X = new float[frame_width * frame_height];
	float* grad_Y = new float[frame_width * frame_height];
	float* grad_T = new float[frame_width * frame_height];

	//初始角点检测
	capture >> frame;
	cv::resize(frame, img, frame_size);
	RGB2GRAY(img, img_GRAY);

	capture >> frame_next;
	cv::resize(frame_next, img_next, frame_size);
	RGB2GRAY(img_next, img_GRAY_next);

	ScharrConvolution(img_GRAY, grad_X, grad_Y);
	getR(result_img, grad_X, grad_Y);
	NMS(result_img, corner_list, 11);
	getBestCorner(corner_list, 20);
	getGradT(img_GRAY, img_GRAY_next, grad_T);
	getUV(grad_X, grad_Y, grad_T, corner_list, velocity_list);
	/*for (int i = 0; i < corner_list.size(); i++) {  //输出各角点速度
		cout << "U:" << velocity_list[i][0] << "V:" << velocity_list[i][1] << endl;
	}*/

	//定义窗口
	mask = Mat::zeros(frame_height, frame_width, img.type());
	namedWindow("img_next", WINDOW_AUTOSIZE);  //显示光流和图像结合的窗口
	imshow("img_next", img_next);
	namedWindow("mask", WINDOW_AUTOSIZE);  //只显示光流的窗口
	imshow("mask", mask);
	namedWindow("origin img", WINDOW_AUTOSIZE);
	imshow("origin img", img);
	

	//角点计算完毕,进入循环
	while (true) {
		//取下一帧图片
		img_GRAY = img_GRAY_next.clone();
		capture >> frame_next;
		if (frame_next.empty()) {
			cout << "视频处理完毕,退出循环" << endl;
			break;
		}
		cv::resize(frame_next, img_next, frame_size);
		imshow("origin img", img_next);

		//在mask上画出光流并根据速度更新坐标
		for (int i = 0; i < corner_list.size(); i++) {
			int y = corner_list[i][0] + 0.5;
			int x = corner_list[i][1] + 0.5;
			int y_next = corner_list[i][0] + velocity_list[i][0] + 0.5;
			int x_next = corner_list[i][1] + velocity_list[i][1] + 0.5;
			circle(img_next, Point(x_next, y_next), 3, Scalar(173, 255, 61));
			line(mask, Point(x, y), Point(x_next, y_next), Scalar(255, 255, 255), 2);
			cout << x << "  " << y << endl;
			corner_list[i][0] += 2 * velocity_list[i][0];
			corner_list[i][1] += 2 * velocity_list[i][1];
		}

		//计算新帧的速度
		RGB2GRAY(img_next, img_GRAY_next);
		ScharrConvolution(img_GRAY_next, grad_X, grad_Y);
		getR(result_img, grad_X, grad_Y);
		getGradT(img_GRAY, img_GRAY_next, grad_T);
		getUV(grad_X, grad_Y, grad_T, corner_list, velocity_list, 10);

		//合并mask和img_next到img_next,显示最终图像
		cv::add(img_next, mask, img_next);
		imshow("mask", mask);
		imshow("img_next", img_next);

		//当点击Esc键时推出循环
		char c = (char)waitKey(10);
		if (c == 27)
			break;
		else if (c == 'p')
			waitKey(0);
	}
	
	//释放动态数组
	delete[] result_img;
	delete[] grad_X;
	delete[] grad_Y;
	delete[] grad_T;

	return 0;
}


/*将彩色图转为灰度图
@param img 待转彩色图
@param GRAY_img 存放灰度图的Mat类型变量,无需初始化
*/
void RGB2GRAY(Mat& img, Mat& GRAY_img) {
	GRAY_img = Mat::zeros(img.size(), CV_8UC1);
	uchar* p_GRAY_img = GRAY_img.data;
	uchar* p_img = img.data;
	for (int i = 0; i < frame_height; i++) {
		for (int j = 0; j < frame_width; j++) {
			p_GRAY_img[i * frame_width + j] = 0.114 * p_img[(i * frame_width + j) * frame_channels] + 0.587 * p_img[(i * frame_width + j) * frame_channels + 1] + 0.299 * p_img[(i * frame_width + j) * frame_channels + 2];
		}
	}
	return;
}


/*使用Scharr算子求XY方向梯度,就是一个卷积的过程,也可以选其他的算子,但注意必须是可以分为两部分分别计算X和Y方向的算子,如Sobel算子、Roberts算子、Prewitt算子
@param img 待求梯度的图像,注意必须是灰度图像
@param grad_X 存放X向梯度的数组指针,数组必须在传入前初始化
@param grad_Y 存放Y向梯度的数组指针,数组必须在传入前初始化
*/
void ScharrConvolution(Mat& img, float* grad_X, float* grad_Y) {
	int ScharrX_operator[3][3] = { -3,0,3,-10,0,10,-3,0,3 };  //Scharr算子
	int ScharrY_operator[3][3] = { -3,-10,3,0,0,0,3,10,3 };
	uchar* p_img = img.data;
	for (int i = 0; i < frame_height - 2; i++) {  //卷积
		for (int j = 0; j < frame_width - 2; j++) {
			float conv_X = 0;
			float conv_Y = 0;
			for (int x = 0; x < 3; x++) {
				for (int y = 0; y < 3; y++) {
					conv_X += p_img[(i + x) * frame_width + j + y] * ScharrX_operator[x][y];
					conv_Y += p_img[(i + x) * frame_width + j + y] * ScharrY_operator[x][y];
				}
			}
			grad_X[(i + 1) * frame_width + j + 1] = abs(conv_X) / 16;
			grad_Y[(i + 1) * frame_width + j + 1] = abs(conv_Y) / 16;
		}
	}
	return;
}


/*获得角点响应值R,用于获取最优角点,其中R=Det-0.05Trace^2
@param result_img 存放各点的角点响应值的数组,需在传入前初始化
@param grad_X 存放X向梯度的数组指针,由ScharrConvolution函数计算
@param grad_Y 存放Y向梯度的数组指针,由ScharrConvolution函数计算
*/
void getR(double* result_img, float* grad_X, float* grad_Y) {
	double gauss[3][3] = { 0.0571, 0.1248, 0.0571, 0.1248, 0.2725, 0.1248, 0.0571, 0.1248, 0.0571 };
	double trace = 0, det = 0;
	double* gradXX = new double[frame_width * frame_height];
	double* gradYY = new double[frame_width * frame_height];
	double* gradXY = new double[frame_width * frame_height];
	for (int i = 1; i < frame_height - 1; i++) {
		for (int j = 1; j < frame_width - 1; j++) {
			gradXX[i * frame_width + j] = grad_X[i * frame_width + j] * grad_X[i * frame_width + j];
			gradYY[i * frame_width + j] = grad_Y[i * frame_width + j] * grad_Y[i * frame_width + j];
			gradXY[i * frame_width + j] = grad_X[i * frame_width + j] * grad_Y[i * frame_width + j];
		}
	}
	for (int i = 1; i < frame_height - 1; i++) {
		for (int j = 1; j < frame_width - 1; j++) {
			double gaussXX = 0, gaussYY = 0, gaussXY = 0;
			for (int x = -1; x < 2; x++) {
				for (int y = -1; y < 2; y++) {
					gaussXX += gauss[x + 1][y + 1] * gradXX[(i + x) * frame_width + j + y];
					gaussYY += gauss[x + 1][y + 1] * gradYY[(i + x) * frame_width + j + y];
					gaussXY += gauss[x + 1][y + 1] * gradXY[(i + x) * frame_width + j + y];
				}
			}
			trace = gaussXX + gaussYY;
			det = gaussXX * gaussYY - gaussXY * gaussXY;
			result_img[i * frame_width + j] = det - 0.05 * trace * trace;
		}
	}
	delete[] gradXX;
	delete[] gradXY;
	delete[] gradYY;
	return;
}


/*非极大值抑制,默认窗口宽度为40
@param result_img 存放各点的角点响应值的数组
@param corner_list 存放角点坐标及对应R值的vector容器,无需初始化
@param core 进行非极大值抑制时的窗口大小,默认为40
*/
void NMS(double* result_img, vector<vector<double>>& corner_list, int core) {
	int semi_core = core / 2;
	for (int i = semi_core + 1; i < frame_height - semi_core; i++) {
		for (int j = semi_core + 1; j < frame_width - semi_core; j++) {
			if (result_img[i * frame_width + j] > 20000) {
				bool is_max = 1;
				// double resulttt = result_img[i * frame_width + j];
				for (int x = -semi_core; x < semi_core + 1; x++) {
					for (int y = -semi_core; y < semi_core + 1; y++) {
						if (result_img[(i + x) * frame_width + j + y] > result_img[i * frame_width + j]) {
							is_max = 0;
						}
					}
				}
				if (is_max) {
					/*cout << i << "," << j << endl;*/
					vector<double> tmp{ double(i), double(j), result_img[i * frame_width + j] };
					corner_list.push_back(tmp);
				}
			}
		}
	}
}


/*自定义排序函数,比较两vector容器中第三个元素(响应值R),用于调用vector排序函数根据R进行排序
@param corner_list1 进行比较的第一个vector容器
@param corner_list2 进行比较的第二个vector容器
*/
bool myCompare(vector<double>& corner_list1, vector<double>& corner_list2) {
	return corner_list1[2] > corner_list2[2];
}


/*根据R值由大到小对坐标进行排序,获取前n个R值最大的坐标,若整体不足或等于n个则全部返回
@param corner_list 存放坐标和对应R值的vector容器
@param n 期望返回的坐标个数
*/
void getBestCorner(vector<vector<double>>& corner_list, int n) {
	if (corner_list.size() <= n) {  //若角点数量不足或等于n,则退出函数
		return;
	}
	sort(corner_list.begin(), corner_list.end(), myCompare);
	corner_list.erase(corner_list.begin() + n, corner_list.end());
	//for (int i = 0; i < corner_list.size(); i++) {  //查看坐标,可在选中后,先后按ctrl+K和ctrl+U来取消注释(Visual Studio)
	//	cout << "x:" << corner_list[i][0] << ",y:" << corner_list[i][1] << ",R:" << corner_list[i][2] << endl;
	//}
	cout << "已提取前" << corner_list.size() << "个角点,目标个数为:" << n << endl;

	return;
}


/*计算t方向梯度
@param img 原图像,需为灰度图
@param img_next 下一帧图像,需为灰度图
@param grad_T 存放T向梯度(两帧差值)的数组指针,需在传入前初始化
*/
void getGradT(Mat img, Mat img_next, float* grad_T) {
	uchar* p_img = img.data;
	uchar* p_img_next = img_next.data;
	for (int i = 0; i < frame_height; i++) {//卷积
		for (int j = 0; j < frame_width; j++) {
			grad_T[i * frame_width + j] = abs(p_img_next[i * frame_width + j] - p_img[i * frame_width + j]);
		}
	}
}


/*计算UV速度:已知 [[XX,XY],[XY,YY]] 矩阵的逆为 1/(XX*YY-XY*XY)[[YY,-XY],[-XY,XX]],
* 则由等式 [U,V]=[[XX,XY],[XY,YY]]逆*(-[XT,YT]) 可以计算出U和V
@param grad_X/grad_Y/grad_T 各向梯度
@param corner_list 存放坐标和对应R值的vector容器
@param velcity_list 存放坐标对应点的速度的vector容器,无需初始化
@param window 计算速度的窗口大小,默认为7,尽量为奇数
*/
void getUV(float* grad_X, float* grad_Y, float* grad_T, vector<vector<double>>& corner_list, vector<vector<double>>& velcity_list, int window) {
	int row, col;  //保存从容器中读取的角点的行和列,对应坐标则为y和x
	velcity_list.clear();  //清除上次迭代计算的速度值
	for (int i = 0;;) {
		row = corner_list[i][0] + 0.5;
		col = corner_list[i][1] + 0.5;
		if (row + window / 2 > frame_height - 1 || row - window / 2 < 0 || col + window / 2 > frame_width - 1 || col - window / 2 < 0) {
			corner_list.erase(corner_list.begin() + i);  //如果清除了第i个元素,则后面的元素会递补第i个位置,所以无需i++
		}
		else {
			i++;
		}
		if (i >= corner_list.size()) {
			break;
		}
	}
	for (int i = 0; i < corner_list.size(); i++) {
		row = corner_list[i][0] + 0.5;
		col = corner_list[i][1] + 0.5;
		double gradXX(0), gradXY(0), gradYY(0), gradXT(0), gradYT(0);

		for (int m = -window / 2; m < window / 2 + 1; m++) {
			for (int n = -window / 2; n < window / 2 + 1; n++) {
				gradXX += grad_X[(row + m) * frame_width + col + n] * grad_X[(row + m) * frame_width + col + n];
				gradXY += grad_X[(row + m) * frame_width + col + n] * grad_Y[(row + m) * frame_width + col + n];
				gradYY += grad_Y[(row + m) * frame_width + col + n] * grad_Y[(row + m) * frame_width + col + n];
				gradXT += grad_X[(row + m) * frame_width + col + n] * grad_T[(row + m) * frame_width + col + n];
				gradYT += grad_Y[(row + m) * frame_width + col + n] * grad_T[(row + m) * frame_width + col + n];
			}
		}
		double deno = gradXX * gradYY - gradXY * gradXY;  //求矩阵逆生成的分母
		double U = (gradYY * gradXT + (-gradXY) * gradYT) / deno;
		double V = ((-gradXY) * gradXT + gradXX * gradYT) / deno;
		vector<double> UV{ U, V };
		velcity_list.push_back(UV);
	}
}

  • 1
    点赞
  • 14
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值