代码 光流法 图像处理

2016-03-18 16:04:36 chentravelling 阅读数 10738

openCV光流法追踪运动物体

email:chentravelling@163.com

一、光流简介

摘自:zouxy09

        光流的概念是Gibson在1950年首先提出来的。它是空间运动物体在观察成像平面上的像素运动的瞬时速度,是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系,从而计算出相邻帧之间物体的运动信息的一种方法。一般而言,光流是由于场景中前景目标本身的移动、相机的运动,或者两者的共同运动所产生的。

研究光流场的目的就是为了从图片序列中近似得到不能直接得到的运动场。运动场,其实就是物体在三维真实世界中的运动;光流场,是运动场在二维图像平面上(人的眼睛或者摄像头)的投影。

       那通俗的讲就是通过一个图片序列,把每张图像中每个像素的运动速度和运动方向找出来就是光流场。那怎么找呢?咱们直观理解肯定是:第t帧的时候A点的位置是(x1, y1),那么我们在第t+1帧的时候再找到A点,假如它的位置是(x2,y2),那么我们就可以确定A点的运动了:(ux, vy) = (x2, y2) - (x1,y1)。

        那怎么知道第t+1帧的时候A点的位置呢? 这就存在很多的光流计算方法了。

       1981年,Horn和Schunck创造性地将二维速度场与灰度相联系,引入光流约束方程,得到光流计算的基本算法。人们基于不同的理论基础提出各种光流计算方法,算法性能各有不同。Barron等人对多种光流计算技术进行了总结,按照理论基础与数学方法的区别把它们分成四种:基于梯度的方法、基于匹配的方法、基于能量的方法、基于相位的方法。近年来神经动力学方法也颇受学者重视。

OpenCV中实现了不少的光流算法。

可参考:http://www.opencv.org.cn/opencvdoc/2.3.2/html/modules/video/doc/motion_analysis_and_object_tracking.html

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光流算法,对运动进行跟踪。

二、代码

#include <opencv2\opencv.hpp>
#include <iostream>

using namespace std;
using namespace cv;
const int MAX_CORNERS = 100;

int main()
{
	IplImage* preImage = cvLoadImage("image133.pgm",CV_LOAD_IMAGE_GRAYSCALE);
	IplImage* curImage = cvLoadImage("image134.pgm",CV_LOAD_IMAGE_GRAYSCALE);
	IplImage* curImage_ = cvLoadImage("image134.pgm");
	CvPoint2D32f* preFeatures = new CvPoint2D32f[ MAX_CORNERS ];//前一帧中特征点坐标(通过cvGoodFeaturesToTrack获得)  
	CvPoint2D32f* curFeatures = new CvPoint2D32f[ MAX_CORNERS ];//在当前帧中的特征点坐标(通过光流法获得)
	double qlevel; //特征检测的指标  
	double minDist;//特征点之间最小容忍距离  
	vector<uchar> status; //特征点被成功跟踪的标志  
	vector<float> err; //跟踪时的特征点小区域误差和 
	CvSize      img_sz    = cvGetSize(preImage);  
	IplImage* eig_image = cvCreateImage( img_sz, IPL_DEPTH_32F, 1 );//缓冲区  
	IplImage* tmp_image = cvCreateImage( img_sz, IPL_DEPTH_32F, 1 );
	int num = MAX_CORNERS;
	cvGoodFeaturesToTrack(//获取特征点
		preImage,
		eig_image,
		tmp_image,
		preFeatures,
		&num,
		0.01,  
		5.0,  
		0,  
		3,  
		0,  
		0.04  
		);
	CvSize pyr_sz = cvSize( curImage->width+8, curImage->height/3 );  
	IplImage* pyrA = cvCreateImage( pyr_sz, IPL_DEPTH_32F, 1 );  
	IplImage* pyrB = cvCreateImage( pyr_sz, IPL_DEPTH_32F, 1 );  
	
	char features_found[ MAX_CORNERS ];  
	float feature_errors[ MAX_CORNERS ];  

	
	cvCalcOpticalFlowPyrLK(//计算光流
		preImage,
		curImage,
		pyrA,
		pyrB,
		preFeatures,
		curFeatures,
		MAX_CORNERS,
		cvSize(10,10),
		5,  
		features_found,  
		feature_errors,  
		cvTermCriteria( CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, .3 ),  
		0  
		);

	
	for(int i=0;i<MAX_CORNERS;i++)//画线
	{
	cvLine(
		curImage_,
		Point(preFeatures[i].x,preFeatures[i].y),
		Point(curFeatures[i].x,curFeatures[i].y),
		CV_RGB(255,0,0),
		4,
		CV_AA,
		0);
	}
	
	cvShowImage("outPutImage",curImage_);
	//cvSaveImage("outPutImage.pgm",curImage);
	waitKey(0);
	return 0;
}

最终结果:



2019-12-26 20:40:11 fhcfhc1112 阅读数 67

1.说文解字
光流法(optical flow)是图像处理/计算机视觉里面一个研究子问题。首先解释一下它的字面意思:optical这个词的中文意思是视觉的,光学的(形容词),flow这个词的中文意思是流(名词),两个单词结合起来被翻译成了光流。但是从理解上我更倾向于把optical 与视觉的联系起来。
生物学上的意义:运动物体在视网膜上的形成的一种感知。这种感知是由物体运动造成的,然后使用我们的语言来描述这种感知——流。
数学物理学上的意义:光流算法的输入是连续的图像序列,输出是运动场。一般的理解说到这就完事了,但是这样的解释依然存在一个阻碍我们理解的词——场。这个词在很多地方都见过,如磁场,量子场,这一个个的甩出来,炸得我们头昏眼花,可谓是相当的神秘。其它场合中场的意思我也不懂,但在这篇文章里面,场是一个函数。那么既然是函数就得有自变量和因变量。自变量是空间位置,因变量是一个二维向量,这个二维向量就表示这个空间位置的速度——速度的大小和方向。这下我们理解了运动场,就可以清晰的理解光流是什么了。
哲学上的意义:那么光流有啥用呢?这就需要我们伟大的马克思主义登场了——世界是物质的,物质是运动的,运动是有规律的,规律是可以认识的,物质的运动包含了规律,包含了信息,于是一切的一切就回到了图像处理的雄心壮志了,我们的征途是通过图像来表示,理解和决策世界。

官方定义:光流是空间运动物体在观察成像平面上的像素运动的瞬时速度(矢量),是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系,从而计算出相邻帧之间像素的运动信息(抽象成运动场也行)的一种方法。

  1. 方法
    相信各位看到这里,对光流就有了一定的认识了。光流法最终要求的目标就是运动场。那么如何求呢?在此介绍两种经典的方法:LK算法和HS算法。
    LK算法假设:
    (1)像素的亮度不变;
    (2)像素的移动缓慢;
    (3)任一像素的附近像素运动信息保持抑制。
2013-09-05 13:41:27 pi9nc 阅读数 1102

2010-02-27 13:50:23|  分类: 默认分类|字号 订阅

光流法的介绍(含C++代码)
2009-07-18 12:31
光流法是比较经典的运动估计方法,本文不仅叙述简单明了,而且附代码,故收藏.
在空间中,运动可以用运动场描述。而在一个图像平面上,物体的运动往往是通过图像序列中不同图象灰度分布的不同体现的。从而,空间中的运动场转移到图像上就表示为光流场,光流场反映了图像上每一点灰度的变化趋势。
光流可以看作带有灰度的像素点在图像平面运动产生的瞬时速度场。下面我们推导光流方程:
假设E(x,y,t)为(x,y)点在时刻t的灰度(照度)。设t+dt时刻该点运动到(x+dx,y+dy)点,他的照度为E(x+dx,y+dy,t+dt)。我们认为,由于对应同一个点,所以
E(x,y,t) = E(x+dx,y+dy,t+dt) —— 光流约束方程
将上式右边做泰勒展开,并令dt->0,则得到:Exu+Eyv+Et = 0,其中:
Ex = dE/dx Ey = dE/dy Et = dE/dt u = dx/dt v = dy/dt
上面的Ex,Ey,Et的计算都很简单,用离散的差分代替导数就可以了。光流法的主要任务就是通过求解光流约束方程求出u,v。但是由于只有一个方程,所以这是个病态问题。所以人们提出了各种其他的约束方程以联立求解。但是由于我们用于摄像机固定的这一特定情况,所以问题可以大大简化。
摄像机固定的情形
在摄像机固定的情形下,运动物体的检测其实就是分离前景和背景的问题。我们知道对于背景,理想情况下,其光流应当为0,只有前景才有光流。所以我们并不要求通过求解光流约束方程求出u,v。我么只要求出亮度梯度方向的速率就可以了,即求出sqrt(u*u+v*v)。
而由光流约束方程可以很容易求到梯度方向的光流速率为 V = abs(Et/sqrt(Ex*Ex+Ey*Ey))。这样我们设定一个阈值T。
V(x,y) > T 则(x,y)是前景 ,反之是背景
C++实现
在实现摄像机固定情况的光流法时,需要有两帧连续的图像,下面的算法针对RGB24格式的图像计算光流:
void calculate(unsigned char* buf)
{
   int Ex,Ey,Et;
   int gray1,gray2;
   int u;
   int i,j;
   memset(opticalflow,0,width*height*sizeof(int));
   memset(output,255,size);
   for(i=2;i<height-2;i++){
for(j=2;j<width-2;j++){
gray1 = int(((int)(buf[(i*width+j)*3])
    +(int)(buf[(i*width+j)*3+1])
    +(int)(buf[(i*width+j)*3+2]))*1.0/3);
gray2 = int(((int)(prevframe[(i*width+j)*3])
    +(int)(prevframe[(i*width+j)*3+1])
    +(int)(prevframe[(i*width+j)*3+2]))*1.0/3);
Et = gray1 - gray2;
gray2 = int(((int)(buf[(i*width+j+1)*3])
    +(int)(buf[(i*width+j+1)*3+1])
    +(int)(buf[(i*width+j+1)*3+2]))*1.0/3);
Ex = gray2 - gray1;
gray2 = int(((int)(buf[((i+1)*width+j)*3])
    +(int)(buf[((i+1)*width+j)*3+1])
    +(int)(buf[((i+1)*width+j)*3+2]))*1.0/3);
Ey = gray2 - gray1;
Ex = ((int)(Ex/10))*10;
Ey = ((int)(Ey/10))*10;
Et = ((int)(Et/10))*10;
u = (int)((Et*10.0)/(sqrt((Ex*Ex+Ey*Ey)*1.0))+0.1);
opticalflow[i*width+j] = u;
if(abs(u)>10){
    output[(i*width+j)*3] = 0;
    output[(i*width+j)*3+1] = 0;
    output[(i*width+j)*3+2] = 0;
}
}
   }
   memcpy(prevframe,buf,size);
}
//////////////////////////////////////////////////////////////////////////
/另一个代码
/
/////////////////////////////////////////////////////////////////////////////
WW_RETURN HumanMotion::ImgOpticalFlow(IplImage *pre_grey,IplImage *grey)
/*************************************************
   Function:
   Description:   光流法计算运动速度与方向    
   Date: 2006-6-14
   Author: 
   Input:                      
   Output:       
   Return:       
   Others:       
*************************************************/
{
IplImage *velx = cvCreateImage( cvSize(grey->width ,grey->height),IPL_DEPTH_32F, 1 );
IplImage *vely = cvCreateImage( cvSize(grey->width ,grey->height),IPL_DEPTH_32F, 1 );
velx->origin =   vely->origin = grey->origin;
CvSize winSize = cvSize(5,5);
cvCalcOpticalFlowLK( prev_grey, grey, winSize, velx, vely );

cvAbsDiff( grey,prev_grey, abs_img );
cvThreshold( abs_img, abs_img, 29, 255, CV_THRESH_BINARY); 
CvScalar xc,yc; 
for(int y =0 ;y<velx->height; y++)
   for(int x =0;x<velx->width;x++ )
   {
xc = cvGetAt(velx,y,x);
yc = cvGetAt(vely,y,x);

float x_shift= (float)xc.val[0]; 
float y_shift= (float)yc.val[0]; 
const int winsize=5;   //计算光流的窗口大小

if((x%(winsize*2)==0) && (y%(winsize*2)==0) ) 
{
if(x_shift!=0 || y_shift!=0)
{

    if(x>winsize && y>winsize && x <(velx->width-winsize) && y<(velx->height-winsize) )
    {
   cvSetImageROI( velx, cvRect( x-winsize, y-winsize, 2*winsize, 2*winsize));
   CvScalar total_x = cvSum(velx);
   float xx = (float)total_x.val[0];
   cvResetImageROI(velx);
   cvSetImageROI( vely, cvRect( x-winsize, y-winsize, 2*winsize, 2*winsize));
   CvScalar total_y = cvSum(vely);
   float yy = (float)total_y.val[0];
   cvResetImageROI(vely);
   
   cvSetImageROI( abs_img, cvRect( x-winsize, y-winsize, 2*winsize, 2*winsize));
   CvScalar total_speed = cvSum(abs_img);
   float ss = (float)total_speed.val[0]/(4*winsize*winsize)/255;
   cvResetImageROI(abs_img);
   const double ZERO = 0.000001;
   const double pi = 3.1415926;
   double alpha_angle;
   if(xx<ZERO && xx>-ZERO)
   alpha_angle = pi/2;
   else
   alpha_angle = abs(atan(yy/xx));
   
   if(xx<0 && yy>0) alpha_angle = pi - alpha_angle ;
   if(xx<0 && yy<0) alpha_angle = pi + alpha_angle ;
   if(xx>0 && yy<0) alpha_angle = 2*pi - alpha_angle ;

   
   CvScalar line_color;
   float scale_factor = ss*100;
   line_color = CV_RGB(255,0,0);
   CvPoint pt1,pt2;
   pt1.x = x; 
   pt1.y = y;
   pt2.x = static_cast<int>(x + scale_factor*cos(alpha_angle));
   pt2.y = static_cast<int>(y + scale_factor*sin(alpha_angle));
   cvLine( image, pt1, pt2 , line_color, 1, CV_AA, 0 );
   CvPoint p;
   p.x = (int) (pt2.x + 6 * cos(alpha_angle - pi / 4*3));
   p.y = (int) (pt2.y + 6 * sin(alpha_angle - pi / 4*3));
   cvLine( image, p, pt2, line_color, 1, CV_AA, 0 );
   p.x = (int) (pt2.x + 6 * cos(alpha_angle + pi / 4*3));
   p.y = (int) (pt2.y + 6 * sin(alpha_angle + pi / 4*3));
   cvLine( image, p, pt2, line_color, 1, CV_AA, 0 );
   /*
   line_color = CV_RGB(255,255,0);
   pt1.x = x-winsize;
   pt1.y = y-winsize;
   pt2.x = x+winsize;
   pt2.y = y+winsize;
   cvRectangle(image, pt1,pt2,line_color,1,CV_AA,0);
   */
    }
}
}
   }

cvShowImage( "Contour", abs_img);
cvShowImage( "Contour2", vely);
cvReleaseImage(&velx);
cvReleaseImage(&vely);
cvWaitKey(20);

return WW_OK;
2019-02-27 18:00:45 weixin_37552816 阅读数 3568

光流法的定义

光流法是空间运动物体在观察成像平面上的像素运动的瞬时速度,是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系,从而计算出相邻帧之间物体的运动信息的一种方法。一般而言,光流是由于场景中前景目标本身的移动、相机的运动,或者两者的共同运动所产生的。
简单来说,光流是空间运动物体在观测成像平面上的像素运动的“瞬时速度”。光流的研究是利用图像序列中的像素强度数据的时域变化和相关性来确定各自像素位置的“运动”。研究光流场的目的就是为了从图片序列中近似得到不能直接得到的运动场。

  光流法的前提假设:
 (1)相邻帧之间的颜色恒定,对于灰度图来说,亮度恒定;
 (2)相邻视频帧的取帧时间连续,或者,相邻帧之间物体的运动比较“微小”;
 (3)保持空间一致性;即,同一子图像的像素点具有相同的运动

上面的三个假设对光流法的原理公式推导很重要,具体的推导公式不在具体阐述。同时也是判断使用光流法的使用条件,光流法主要分为两种,计算部分像素运动的:稀疏光流,以Lucas-Kanade为代表,计算所有像素运动的:稠密光流。下面将结合代码具体理解。

Lucas-Kanade光流原理

L-K算法基于上面三个假设。
(1)亮度恒定,就是同一点随着时间的变化,其亮度不会发生改变。这是基本光流法的假定(所有光流法变种都必须满足),用于得到光流法基本方程;
(2)小运动,这个也必须满足,就是时间的变化不会引起位置的剧烈变化,这样灰度才能对位置求偏导(换句话说,小运动情况下我们才能用前后帧之间单位位置变化引起的灰度变化去近似灰度对位置的偏导数),这也是光流法不可或缺的假定;
(3)空间一致,一个场景上邻近的点投影到图像上也是邻近点,且邻近点速度一致。这是Lucas-Kanade光流法特有的假定,因为光流法基本方程约束只有一个,而要求x,y方向的速度,有两个未知变量。我们假定特征点邻域内做相似运动,就可以连立区域内个个像素点的多个方程求取x,y方向的速度。
这样假设的具体原因可参考:https://blog.csdn.net/zyazky/article/details/52175069

L-K光流法是稀疏光流法,需要找到视频每帧的特征点,这个常用的有Shi-Tomasi角点和SIFT特征点。

Shi-Tomasi角点检测

opencv中调用goodFeaturesToTrack函数来实现角点检测

void goodFeaturesToTrack( InputArray image, OutputArray corners,
                                     int maxCorners, double qualityLevel, double minDistance,
                                     InputArray mask=noArray(), int blockSize=3,
                                     bool useHarrisDetector=false, double k=0.04 );

第一个参数image:8位或32位单通道灰度图像;
第二个参数corners:位置点向量,保存的是检测到的角点的坐标;

第三个参数maxCorners:定义可以检测到的角点的数量的最大值;

第四个参数qualityLevel:检测到的角点的质量等级,角点特征值小于qualityLevel*最大特征值的点将被舍弃;

第五个参数minDistance:两个角点间最小间距,以像素为单位;

第六个参数mask:指定检测区域,若检测整幅图像,mask置为空Mat();

第七个参数blockSize:计算协方差矩阵时窗口大小;默认为3

第八个参数useHarrisDetector:是否使用Harris角点检测,为false,则使用Shi-Tomasi算子;

第九个参数k:留给Harris角点检测算子用的中间参数,一般取经验值0.04~0.06。第八个参数为false时,该参数不起作用;

mask可以指定角点产生的区域,划定检测范围

python-opencv代码demo

opencv中提供了一个函数来实现K-L算法:cv2.calcOpticalFlowPyrLK()

nextPts,status,err = cv.calcOpticalFlowPyrLK(   prevImg, nextImg, prevPts, nextPts[, status[, err[, winSize[, maxLevel[, criteria[, flags[, minEigThreshold]]]]]]])

返回值:

nextPtrs:输出一个二维点的向量,这个向量可以是用来作为光流算法的输入特征点,也是光流算法在当前帧找到特征点的新位置(浮点数)
status 标志:在当前帧当中发现的特征点标志status==1,否则为0
err :向量中的每个特征对应的错误率

输入值:
prevImg :上一帧图片
nextImg :当前帧图片
prevPts :上一帧找到的特征点向量
nextPts :与返回值中的nextPtrs相同
status :与返回的status相同
err :与返回的err相同
winSize :在计算局部连续运动的窗口尺寸(在图像金字塔中)
maxLevel: 图像金字塔层数,0表示不使用金字塔
criteria :寻找光流迭代终止的条件
flags :有两个宏,表示两种计算方法,分别是OPTFLOW_USE_INITIAL_FLOW:表示使用估计值作为寻找到的初始光流,OPTFLOW_LK_GET_MIN_EIGENVALS表示:使用最小特征值作为误差测量
minEigThreshold :该算法计算光流方程的2×2规范化矩阵的最小特征值,除以窗口中的像素数; 如果此值小于minEigThreshold,则会过滤掉相应的功能并且不会处理该光流,因此它允许删除坏点并获得性能提升。

下面附上代码:
读取视频第一帧,检测角点,然后使用K-L算法来迭代跟踪每个角点在每一帧的位置信息,最后画出运动轨迹。在做煤矿安全生产的项目时,考虑过使用光流法来检测传送到运动方向和速度,但是视频灰度图像中,运煤的传送带的灰度图亮度差异太小,导致角点产生了跳变。

import numpy as np
import cv2
import matplotlib.pyplot as plt

cap = cv2.VideoCapture('1/12311sc.mp4')
# params for ShiTomasi corner detection 特征点检测
feature_params = dict( maxCorners = 10,
                       qualityLevel = 0.1,
                       minDistance = 10,
                       blockSize = 3 )

# Parameters for lucas kanade optical flow光流法参数
lk_params = dict(winSize = (15,15),
                  maxLevel = 0,
                  criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))

# Create some random colors 画轨迹
color = np.random.randint(0,255,(100,3))

# Take first frame and find corners in it
ret, old_frame = cap.read()
old_gray = cv2.cvtColor(old_frame, cv2.COLOR_BGR2GRAY)
roi = np.zeros_like(old_gray)
x,y,w,h = 266,143,150,150
roi[y:y+h, x:x+w] = 255
p0 = cv2.goodFeaturesToTrack(old_gray, mask = roi, **feature_params)

# Create a mask image for drawing purposes
mask = np.zeros_like(old_frame)

while(1):
    ret,frame = cap.read()
    if not ret:
      break
    frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

    # calculate optical flow
    p1, st, err = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)

    # Select good points
    good_new = p1[st==1]
    good_old = p0[st==1]

    # draw the tracks
    for i,(new,old) in enumerate(zip(good_new,good_old)):
        a,b = new.ravel()
        c,d = old.ravel()
        mask = cv2.line(mask, (a,b),(c,d), color[i].tolist(), 2)
        frame = cv2.circle(frame,(a,b),3,color[i].tolist(),-1)
    img = cv2.add(frame,mask)

    cv2.imshow('frame',img)
    key = cv2.waitKey(60) & 0xff
    if key == 27:  # 按下ESC时,退出
        break
    elif key == ord(' '):  # 按下空格键时,暂停
        cv2.waitKey(0)

    # Now update the previous frame and previous points
    old_gray = frame_gray.copy()
    p0 = good_new.reshape(-1,1,2)

cv2.destroyAllWindows()
cap.release()

如有疑问或错误,欢迎大家指正交流。

2015-09-23 14:17:48 mydear_11000 阅读数 3098
原创文章,转贴请注明:http://blog.csdn.net/crzy_sparrow/article/details/7407604

   

本文目录:

      一.基于特征点的目标跟踪的一般方法

      二.光流法

      三.opencv中的光流法函数

      四.用类封装基于光流法的目标跟踪方法

      五.完整代码

      六.参考文献


一.基于特征点的目标跟踪的一般方法

      基于特征点的跟踪算法大致可以分为两个步骤:

      1)探测当前帧的特征点;

      2)通过当前帧和下一帧灰度比较,估计当前帧特征点在下一帧的位置;

      3)过滤位置不变的特征点,余下的点就是目标了。

      很显然,基于特征点的目标跟踪算法和1),2)两个步骤有关。特征点可以是Harris角点(见我的另外一篇博文),也可以是边缘点等等,而估计下一帧位置的方法也有不少,比如这里要讲的光流法,也可以是卡尔曼滤波法(咱是控制系的,上课经常遇到这个,所以看光流法看着看着就想到这个了)。

      本文中,用改进的Harris角点提取特征点(见我另一篇博文:http://blog.csdn.net/crzy_sparrow/article/details/7391511),用Lucas-Kanade光流法实现目标跟踪。


二.光流法

      这一部分《learing opencv》一书的第10章Lucas-Kanade光流部分写得非常详细,推荐大家看书。我这里也粘帖一些选自书中的内容。

      另外我对这一部分附上一些个人的看法(谬误之处还望不吝指正):

      1.首先是假设条件:

       (1)亮度恒定,就是同一点随着时间的变化,其亮度不会发生改变。这是基本光流法的假定(所有光流法变种都必须满足),用于得到光流法基本方程;

       (2)小运动,这个也必须满足,就是时间的变化不会引起位置的剧烈变化,这样灰度才能对位置求偏导(换句话说,小运动情况下我们才能用前后帧之间单位位置变化引起的灰度变化去近似灰度对位置的偏导数),这也是光流法不可或缺的假定;

       (3)空间一致,一个场景上邻近的点投影到图像上也是邻近点,且邻近点速度一致。这是Lucas-Kanade光流法特有的假定,因为光流法基本方程约束只有一个,而要求x,y方向的速度,有两个未知变量。我们假定特征点邻域内做相似运动,就可以连立n多个方程求取x,y方向的速度(n为特征点邻域总点数,包括该特征点)。

      2.方程求解

      多个方程求两个未知变量,又是线性方程,很容易就想到用最小二乘法,事实上opencv也是这么做的。其中,最小误差平方和为最优化指标。

      3.好吧,前面说到了小运动这个假定,聪明的你肯定很不爽了,目标速度很快那这货不是二掉了。幸运的是多尺度能解决这个问题。首先,对每一帧建立一个高斯金字塔,最大尺度图片在最顶层,原始图片在底层。然后,从顶层开始估计下一帧所在位置,作为下一层的初始位置,沿着金字塔向下搜索,重复估计动作,直到到达金字塔的底层。聪明的你肯定发现了:这样搜索不仅可以解决大运动目标跟踪,也可以一定程度上解决孔径问题(相同大小的窗口能覆盖大尺度图片上尽量多的角点,而这些角点无法在原始图片上被覆盖)。




三.opencv中的光流法函数

      opencv2.3.1中已经实现了基于光流法的特征点位置估计函数(当前帧位置已知,前后帧灰度已知),介绍如下(摘自opencv2.3.1参考手册):

  1. calcOpticalFlowPyrLK  
  2. Calculates an optical flow for a sparse feature set using the iterative Lucas-Kanade method with pyramids.  
  3.   
  4. void calcOpticalFlowPyrLK(InputArray prevImg, InputArray nextImg, InputArray prevPts,  
  5. InputOutputArray nextPts, OutputArray status, OutputArray err,  
  6. Size winSize=Size(15,15), int maxLevel=3, TermCriteria crite-  
  7. ria=TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 30, 0.01),  
  8. double derivLambda=0.5, int flags=0 )  
  9.   
  10. Parameters  
  11. prevImg – First 8-bit single-channel or 3-channel input image.  
  12. nextImg – Second input image of the same size and the same type as prevImg .  
  13. prevPts – Vector of 2D points for which the flow needs to be found. The point coordinates  
  14. must be single-precision floating-point numbers.  
  15. nextPts – Output vector of 2D points (with single-precision floating-point coordinates)  
  16. containing the calculated new positions of input features in the second image. When  
  17. OPTFLOW_USE_INITIAL_FLOW flag is passed, the vector must have the same size as in the  
  18. input.  
  19. status – Output status vector. Each element of the vector is set to 1 if the flow for the  
  20. corresponding features has been found. Otherwise, it is set to 0.  
  21. err – Output vector that contains the difference between patches around the original and  
  22. moved points.  
  23. winSize – Size of the search window at each pyramid level.  
  24.   
  25. maxLevel – 0-based maximal pyramid level number. If set to 0, pyramids are not used  
  26. (single level). If set to 1, two levels are used, and so on.  
  27. criteria – Parameter specifying the termination criteria of the iterative search algorithm  
  28. (after the specified maximum number of iterations criteria.maxCount or when the search  
  29. window moves by less than criteria.epsilon .  
  30. derivLambda – Not used.  
  31.   
  32. flags – Operation flags:  
  33. – OPTFLOW_USE_INITIAL_FLOW Use initial estimations stored in nextPts . If the  
  34. flag is not set, then prevPts is copied to nextPts and is considered as the initial estimate.  

四.用类封装基于光流法的目标跟踪方法

      废话少说,附上代码,包括特征点提取,跟踪特征点,标记特征点等。

  1. <span style="font-size:18px;">//帧处理基类  
  2. class FrameProcessor{  
  3.     public:  
  4.         virtual void process(Mat &input,Mat &ouput)=0;  
  5. };  
  6.   
  7. //特征跟踪类,继承自帧处理基类  
  8. class FeatureTracker :  public FrameProcessor{  
  9.     Mat gray;  //当前灰度图  
  10.     Mat gray_prev;  //之前的灰度图  
  11.     vector<Point2f> points[2];//前后两帧的特征点  
  12.     vector<Point2f> initial;//初始特征点  
  13.     vector<Point2f> features;//检测到的特征  
  14.     int max_count; //要跟踪特征的最大数目  
  15.     double qlevel; //特征检测的指标  
  16.     double minDist;//特征点之间最小容忍距离  
  17.     vector<uchar> status; //特征点被成功跟踪的标志  
  18.     vector<float> err; //跟踪时的特征点小区域误差和  
  19. public:  
  20.     FeatureTracker():max_count(500),qlevel(0.01),minDist(10.){}  
  21.     void process(Mat &frame,Mat &output){  
  22.         //得到灰度图  
  23.         cvtColor (frame,gray,CV_BGR2GRAY);  
  24.         frame.copyTo (output);  
  25.         //特征点太少了,重新检测特征点  
  26.         if(addNewPoint()){  
  27.             detectFeaturePoint ();  
  28.             //插入检测到的特征点  
  29.             points[0].insert (points[0].end (),features.begin (),features.end ());  
  30.             initial.insert (initial.end (),features.begin (),features.end ());  
  31.         }  
  32.         //第一帧  
  33.         if(gray_prev.empty ()){  
  34.                 gray.copyTo (gray_prev);  
  35.         }  
  36.         //根据前后两帧灰度图估计前一帧特征点在当前帧的位置  
  37.         //默认窗口是15*15  
  38.         calcOpticalFlowPyrLK (  
  39.                 gray_prev,//前一帧灰度图  
  40.                 gray,//当前帧灰度图  
  41.                 points[0],//前一帧特征点位置  
  42.                 points[1],//当前帧特征点位置  
  43.                 status,//特征点被成功跟踪的标志  
  44.                 err);//前一帧特征点点小区域和当前特征点小区域间的差,根据差的大小可删除那些运动变化剧烈的点  
  45.         int k = 0;  
  46.         //去除那些未移动的特征点  
  47.         for(int i=0;i<points[1].size ();i++){  
  48.             if(acceptTrackedPoint (i)){  
  49.                 initial[k]=initial[i];  
  50.                 points[1][k++] = points[1][i];  
  51.             }  
  52.         }  
  53.         points[1].resize (k);  
  54.         initial.resize (k);  
  55.         //标记被跟踪的特征点  
  56.         handleTrackedPoint (frame,output);  
  57.         //为下一帧跟踪初始化特征点集和灰度图像  
  58.         std::swap(points[1],points[0]);  
  59.         cv::swap(gray_prev,gray);  
  60.     }  
  61.   
  62.     void detectFeaturePoint(){  
  63.         goodFeaturesToTrack (gray,//输入图片  
  64.                                  features,//输出特征点  
  65.                                  max_count,//特征点最大数目  
  66.                                  qlevel,//质量指标  
  67.                                  minDist);//最小容忍距离  
  68.     }  
  69.     bool addNewPoint(){  
  70.         //若特征点数目少于10,则决定添加特征点  
  71.         return points[0].size ()<=10;  
  72.     }  
  73.   
  74.     //若特征点在前后两帧移动了,则认为该点是目标点,且可被跟踪  
  75.     bool acceptTrackedPoint(int i){  
  76.         return status[i]&&  
  77.                 (abs(points[0][i].x-points[1][i].x)+  
  78.                   abs(points[0][i].y-points[1][i].y) >2);  
  79.     }  
  80.   
  81.     //画特征点  
  82.     void  handleTrackedPoint(Mat &frame,Mat &output){  
  83.             for(int i=0;i<points[i].size ();i++){  
  84.                 //当前特征点到初始位置用直线表示  
  85.                 line(output,initial[i],points[1][i],Scalar::all (0));  
  86.                 //当前位置用圈标出  
  87.                 circle(output,points[1][i],3,Scalar::all(0),(-1));  
  88.             }  
  89.         }  
  90. };</span>  

五.完整代码
  完整的运行代码有300+行,粘上来太多了,大家去我传的资源里下载吧。

  下载地址:http://download.csdn.net/detail/crzy_sparrow/4183674

  运行结果:



六.参考文献

【1】The classic article by B. Lucas and T. Kanade, An iterative image registration technique with an application to stereo vision in Int. Joint Conference in Artificial Intelligence, pp. 674-679,1981, that describes the original feature point tracking algorithm.
【2】The article by J. Shi and C. Tomasi, Good Features to Track in IEEE Conference on Computer Vision and Pattern Recognition, pp. 593-600, 1994, that describes an improved version of the original feature point tracking algorithm.