多帧雷达回波光流追踪

多帧雷达回波光流追踪

1. 原理

LK光流法计算的是图像之间角点之间的光流(简单来讲就是两帧图像之间匹配到的角点之间的移动矢量)，由于两帧之间角点可能较少，因此计算出来的光流较为稀疏，若使用插值将稀疏光流插值成稠密光流，则误差较大。而多帧图像光流则能较好的减少这一误差，同时保留好图像之间的时间信息。

2. 工作流程

1. 首先计算多帧图像的角点。
2. 通过LK光流法计算两帧图像之间角点光流
3. 去聚类算法删除角点较近的光流
4. 通过反距离权重插值(IDW插值)将稀疏的角点光流插值成稠密光流

下图展示了三帧雷达回波图像的光流追踪：

3. 具体编码实现：

由上文可知，主要需要算法包括角点计算、LK光流算法、去聚类算法、IDW算法四个部分。

3.1 角点计算

void cv::goodFeaturesToTrack( InputArray _image, OutputArray _corners,
                  int maxCorners, double qualityLevel, double minDistance,
                   InputArray _mask, int blockSize,
                     bool useHarrisDetector, double harrisK )

具体参数见：OpenCV: Feature Detection，或百度其他文档。

3.2 LK光流算法

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 
)	

具体参数见:OpenCV: Object Tracking，或百度其他文档。

3.3 去聚类算法

去聚类即去除位置较近的点，取这些点的中值即可。

去聚类能方便后续IDW插值，同时可以有效保证高质量光流

编码如下：

void decluster(vector<Point2f> xy,vector<Point2f> uv,vector<Point2f> &res_xy,vector<Point2f> &res_uv,int scale)
{
	vector<Point2f> f_xy=xy;
	//将位置进行聚类分析。
	for(int i=0;i<xy.size();i++)
	{
		f_xy[i].x=floor(xy[i].x/scale);
		f_xy[i].y=floor(xy[i].y/scale);
	}

	//排序且删除重复向量
	vector<Point2f> u_xy=f_xy;
	sort(u_xy.begin(),u_xy.end(),[](Point2f pts1,Point2f pts2){return pts1.x<pts2.x;});
	for(int i=0;i<u_xy.size();i++)
	{
		for(int j=i+1;j<u_xy.size();)
		{
			if ( (u_xy[i].x == u_xy[j].x)&&(u_xy[i].y==u_xy[j].y) )
			{
				u_xy.erase(u_xy.begin()+j);
			}
			else
			{
				j++;
			}
		}
	}

	//在f_xy中找到与u_xy[i]相同的点，在找到的所有点中取中值。记录并返回
	for(int i=0;i<u_xy.size();i++)
	{
		vector<int> idx;
		for(int j=0;j<f_xy.size();j++)
		{
			if(f_xy[j]==u_xy[i])
			{
				idx.push_back(j);
			}
		}
		Point2f sum_coord{0,0};
		Point2f sum_uv{0,0};
		for(int j=0;j<idx.size();j++)
		{
			sum_coord+=xy[idx[j]];
			sum_uv+=uv[idx[j]];
		}
		Point2f median_coord{0,0};
		median_coord.x=floor(sum_coord.x/idx.size());
		median_coord.y=floor(sum_coord.y/idx.size());
		Point2f median_uv{0,0};
		median_uv.x=sum_uv.x/idx.size();
		median_uv.y=sum_uv.y/idx.size();
		res_xy.push_back(median_coord);
		res_uv.push_back(median_uv);
	}
}

函数API接口如下：

参数	参数类型	描述
xy	vector	输入参数，光流的位置信息，储存了光流的x和y坐标
uv	vector	输入参数，光流的速度信息，储存了光流x和y轴的速度信息
res_xy	vector	输出参数，去聚类以后光流的位置信息。
res_uv	vector	输出参数，去聚类以后光流的速度信息
scale	int	输入参数，去聚类的范围。

3.4 IDW插值

原理较为简单，之后会专门写一篇关于图像插值的博客。现阶段如果不懂可以百度。

以下为光流IDW插值的编码：

vector<float> Distance(Point2f p,vector<Point2f> map,float distance_offset)
{
	vector<float> distance;
	for(int i=0;i<map.size();i++)
	{
		distance.push_back(sqrt(pow(map[i].x-p.x,2)+pow(map[i].y-p.y,2))+distance_offset);
	}
	return distance;
}

void IDW_interpolate(vector<Point2f> xy,vector<Point2f> uv,vector<Point2f> &uv_res,int Height,int Width)
{
	vector<float> dense_flow;
	dense_flow.resize(Height*Width);
	vector<vector<float>> AllDistance;//储存像素点之间的距离
	vector<vector<float>> weights;
	for(int i=0;i<Height;i++)
	{
		for(int j=0;j<Width;j++)
		{
			Point2f p;
			p.y=i;p.x=j;
			AllDistance.push_back(Distance(p,xy,0.5));
		}
	}
	for(int i=0;i<Height;i++)
	{
		for(int j=0;j<Width;j++)
		{
			float sumx=0;
			float sumy=0;
			float weights_sum=0;
			int idx=i*Width+j;
			for(int n=0;n<AllDistance[idx].size();n++)	
			{
				sumx += uv[n].x / (pow(AllDistance[idx][n], 2));
				sumy += uv[n].y / (pow(AllDistance[idx][n], 2));
				weights_sum += 1 / (pow(AllDistance[idx][n], 2));
			}
			Point2f p;
			p.x = sumx / weights_sum;
			p.y = sumy / weights_sum;
			uv_res.push_back(p);
		}
	}
}

3.5 最后

多帧雷达回波光流追踪整个过程就是角点追踪计算成稀疏光流，将光流数据处理以后插值成稠密光流的过程。编码如下：

vector<Point2f> dense_LK(vector<Mat> input_images)
{
	int ImgNum=input_images.size();
	//批量处理多张图片
	vector<Point2f> xy;
	vector<Point2f> uv;
	for(int n=0;n<ImgNum-1;n++)
	{
		Mat prevs_img=input_images[n];
		Mat next_img=input_images[n+1];

		//形态学开运算去除小噪声
		Mat element=getStructuringElement(MORPH_RECT,Size(3,3));
		morphologyEx(prevs_img,prevs_img,MORPH_OPEN,element);
		morphologyEx(next_img,next_img,MORPH_OPEN,element);

		//特征点搜索
		Mat points;
		int max_corners=1000;
		double quality_level=0.01;
		double min_distance=10;
		int blockSize=3;
		bool useHarrisDetector=false;
		double k=0.04;
		goodFeaturesToTrack(prevs_img,points,max_corners,quality_level,min_distance);

		//LK光流跟踪
		Size winSize{21,21};
		int maxLevel=3;
		Mat points1;
		Mat err;
		Mat status;
		calcOpticalFlowPyrLK(prevs_img,next_img,points,points1,status,err,winSize,maxLevel);


		//数据转换
		vector<Point2f> vec_points=Mat_Point2f2vector(points);
		vector<Point2f> vec_points1=Mat_Point2f2vector(points1);
		vector<unsigned char> vec_stt=Mat_uchar2vector(status);

		//去除无效追踪点
		for(int i=0;i<vec_stt.size();)
		{
			if(!vec_stt[i])
			{
				vec_points.erase(vec_points.begin()+i);
				vec_points1.erase(vec_points1.begin()+i);
			}
			else
			{
				i++;
			}
		}

		//保存数据
		vector<Point2f> _xy=vec_points;
		vector<Point2f> _uv;
		for(int i=0;i<vec_points.size();i++)
		{
			_uv.push_back(vec_points1[i]-vec_points[i]);
		}
		for(int i=0;i<_xy.size();i++)
		{
			xy.push_back(_xy[i]);
			uv.push_back(_uv[i]);
		}
	}
	//去聚类稀疏光流矢量
	vector<Point2f> uv_decluster;
	vector<Point2f> xy_decluster;
	decluster(xy,uv,xy_decluster,uv_decluster,10);

	//插值成稠密光流
	vector<Point2f> uv_res;
	int Height=input_images[0].rows;
	int Width=input_images[0].cols;
	IDW_interpolate(xy_decluster,uv_decluster,uv_res,input_images[0].rows,input_images[0].cols);
    return uv_res;
}

不足与改进

1. 整个过程其实并没有很好的保留文章开头所说的时间信息，我认为可以在去聚类这里做点文章，比如加上时间权重来保留所需要的光流(根据需求具体改进）。
2. 该算法并没有优于Opencv中的稠密光流算法，仅做学习过程中的记录。

最后给个关注和点赞，谢谢。