OpenCV高版本已经把粒子滤波的CV方面的condensation算法给去掉了，以前学的condensation算法不能用C++开发还是只能用C版本，（OpenCV3学习（12.4） 粒子滤波Condensation算法）。
粒子滤波其实有很多变种，Rob Hess实现的这种应该是最基本的一种，Sampling Importance Resampling (SIR)，根据重要性重采样。
1）初始化阶段-提取跟踪目标特征
该阶段要人工指定跟踪目标，程序计算跟踪目标的特征，比如可以采用目标的颜色特征。具体到Rob Hess的代码，开始时需要人工用鼠标拖动出一个跟踪区域，然后程序自动计算该区域色调(Hue)空间的直方图，即为目标的特征。直方图可以用一个向量来表示，所以目标特征就是一个N*1的向量V。
2）搜索阶段-放狗
好，我们已经掌握了目标的特征，下面放出很多条狗，去搜索目标对象，这里的狗就是粒子particle。狗有很多种放法。比如，a)均匀的放：即在整个图像平面均匀的撒粒子(uniform distribution)；b)在上一帧得到的目标附近按照高斯分布来放，可以理解成，靠近目标的地方多放，远离目标的地方少放。Rob Hess的代码用的是后一种方法。狗放出去后，每条狗怎么搜索目标呢？就是按照初始化阶段得到的目标特征(色调直方图，向量V)。每条狗计算它所处的位置处图像的颜色特征，得到一个色调直方图，向量Vi，计算该直方图与目标直方图的相似性。相似性有多种度量，最简单的一种是计算sum(abs(Vi-V)).每条狗算出相似度后再做一次归一化，使得所有的狗得到的相似度加起来等于1.
3）决策阶段
我们放出去的一条条聪明的狗向我们发回报告，“一号狗处图像与目标的相似度是0.3”,“二号狗处图像与目标的相似度是0.02”,“三号狗处图像与目标的相似度是0.0003”,“N号狗处图像与目标的相似度是0.013”...那么目标究竟最可能在哪里呢？我们做次加权平均吧。设N号狗的图像像素坐标是(Xn,Yn),它报告的相似度是Wn,于是目标最可能的像素坐标X = sum(Xn*Wn),Y = sum(Yn*Wn).
4）重采样阶段Resampling
既然我们是在做目标跟踪，一般说来，目标是跑来跑去乱动的。在新的一帧图像里，目标可能在哪里呢？还是让我们放狗搜索吧。但现在应该怎样放狗呢？让我们重温下狗狗们的报告吧。“一号狗处图像与目标的相似度是0.3”,“二号狗处图像与目标的相似度是0.02”,“三号狗处图像与目标的相似度是0.0003”,“N号狗处图像与目标的相似度是0.013”...综合所有狗的报告，一号狗处的相似度最高，三号狗处的相似度最低，于是我们要重新分布警力，正所谓好钢用在刀刃上，我们在相似度最高的狗那里放更多条狗，在相似度最低的狗那里少放狗，甚至把原来那条狗也撤回来。这就是Sampling Importance Resampling，根据重要性重采样(更具重要性重新放狗)。
(2)->(3)->(4)->(2)如是反复循环，即完成了目标的动态跟踪。
粒子滤波的核心思想是随机采样+重要性重采样。既然我不知道目标在哪里，那我就随机的撒粒子吧。撒完粒子后，根据特征相似度计算每个粒子的重要性，然后在重要的地方多撒粒子，不重要的地方少撒粒子。所以说粒子滤波较之蒙特卡洛滤波，计算量较小。这个思想和RANSAC算法真是不谋而合。RANSAC的思想也是(比如用在最简单的直线拟合上)，既然我不知道直线方程是什么，那我就随机的取两个点先算个直线出来，然后再看有多少点符合我的这条直线。哪条直线能获得最多的点的支持，哪条直线就是目标直线。想法非常简单，但效果很好。
from:https://blog.csdn.net/yangtrees/article/details/8078612
转载自： http://www.cnblogs.com/tornadomeet/archive/2012/06/23/2559193.html
Opencv实现粒子滤波算法
目标跟踪过程分为2部分，即目标特征提取和目标跟踪算法。
其中目标特征提取又包括以下几种：1. 各种色彩空间直方图，利用色彩空间的直方图分布作为目标跟踪的特征，可以减少物体远近距离的影响，因为其颜色分布大致相同。2.轮廓特征，提取目标的轮廓特征，可以加快算法的速度，且可以在目标有小部分影响的情况下同样有效果。3. 纹理特征，如果被跟踪目标是有纹理的，则根据其纹理特征来跟踪效果会有所改善。
目标跟踪算法目前大概分为以下4种：1. 基于meanshift算法，即利用meanshift算法可以快速找到领域目标最相似的地方，效果还不错，但是其只能找到局部最大值，且不能解决遮挡问题以及不能自适应跟踪目标的形状，方向等。其后面有学者对其做了改进，比如说camshift，就可以自适应物体的大小，方向，具有较好的跟踪效果。2. Kalman滤波的思想，该思想是利用物体的运动模型来，即服从高斯模型，来对目标状态进行预测，然后与观察模型进行比较，根据2者之间的误差来寻找运动目标的状态，但是该算法的精度不高，因为其高斯运动模型在现实生活中很多条件下并得不到满足，并且该算法对杂乱的背景也很敏感。3. 基于粒子滤波的思想，每次通过实验可以重采样粒子的分布，根据该分布对粒子进行扩散，然后通过扩散的结果来观察目标的状态，最后更新目标的状态。该算法最大的特点是跟踪速度快，且能解决一部分遮挡问题，在实际应用过程中越来越多。4.基于目标建模的方法。该方法具有一定的针对性，需要提前知道所需跟踪的目标是什么，比如说车辆，人脸，行人等。由于已经知道了跟踪目标，所以必须对目标进行建模，然后利用该模型来进行跟踪。该方法的局限性是必须提前知道所跟踪的目标是什么，因而其推广性比较差。
实现过程
1. 在摄像头采集到的视频序列中手动标注一个目标区域，用矩形框表示，并计算该目标区域内的直方图，作为匹配模板。
2. 在选中的目标区域中心处撒预定值为100的粒子数目，并对这100个粒子的结构体都初始化为同样的值，比如说粒子位置为目标区域中心，粒子矩形长宽为目标矩形长宽，粒子的初始尺寸为1.等等。
其粒子的结构体和注释如下：
/****定义粒子结构体****/
typedef struct particle
{
int orix,oriy;//原始粒子坐标
int x,y;//当前粒子的坐标
double scale;//当前粒子窗口的尺寸
int prex,prey;//上一帧粒子的坐标
double prescale;//上一帧粒子窗口的尺寸
Rect rect;//当前粒子矩形窗口
Mat hist;//当前粒子窗口直方图特征
double weight;//当前粒子权值
}PARTICLE;3. 利用二阶动态模型对这100个粒子进行随机扩散，每个粒子根据二阶随机动态模型扩散到一个新的位置。
4. 计算以新位置处每个粒子为中心的矩形框内图像的直方图特征，然后每个特征都与目标模板直方图进行比较，计算其相似度。
5. 根据计算得到的相似度算出每个粒子的权值，即相似度伟大的权值越大，并且对权值进行归一化。
6. 取最大权值处的粒子中心为跟踪目标中心，并根据粒子结构体中的尺寸成员算出目标的矩形框。
7. 根据上一个状态的粒子分布情况，按照其权值乘以粒子总数100进行重采用，即权值大的粒子其采样得到的粒子数也大。这样可以重采用重新得到100个粒子。
8. 根据这100个粒子的情况，重新进入步骤3进行循环，经过对粒子权值分布的重采样，根据二阶随机运动模型进行粒子扩散，根据直方图的巴氏距离推断粒子权值，然后取权值最大的那个作为目标，如此循环。
实例：
// particle_tracking.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include <opencv2/core/core.hpp>
#include "opencv2/imgproc/imgproc.hpp"
#include <opencv2/highgui/highgui.hpp>
#include <stdio.h>
#include <iostream>
using namespace cv;
using namespace std;
Rect select;
bool select_flag=false;
bool tracking=false;//跟踪标志位
bool select_show=false;
Point origin;
Mat frame,hsv;
int after_select_frames=0;//选择矩形区域完后的帧计数
/****rgb空间用到的变量****/
//int hist_size[]={16,16,16};//rgb空间各维度的bin个数
//float rrange[]={0,255.0};
//float grange[]={0,255.0};
//float brange[]={0,255.0};
//const float *ranges[] ={rrange,grange,brange};//range相当于一个二维数组指针
/****hsv空间用到的变量****/
int hist_size[]={16,16,16};
float hrange[]={0,180.0};
float srange[]={0,256.0};
float vrange[]={0,256.0};
//int hist_size[]={32,32,32};
//float hrange[]={0,359.0.0};
//float srange[]={0,1.0};
//float vrange[]={0,1.0};
const float *ranges[]={hrange,srange,vrange};
int channels[]={0,1,2};
/****有关粒子窗口变化用到的相关变量****/
int A1=2;
int A2=-1;
int B0=1;
double sigmax=1.0;
double sigmay=0.5;
double sigmas=0.001;
/****定义使用粒子数目宏****/
#define PARTICLE_NUMBER 100 //如果这个数设定太大，经测试这个数字超过25就会报错，则在运行时将会出现错误
/****定义粒子结构体****/
typedef struct particle
{
int orix,oriy;//原始粒子坐标
int x,y;//当前粒子的坐标
double scale;//当前粒子窗口的尺寸
int prex,prey;//上一帧粒子的坐标
double prescale;//上一帧粒子窗口的尺寸
Rect rect;//当前粒子矩形窗口
Mat hist;//当前粒子窗口直方图特征
double weight;//当前粒子权值
}PARTICLE;
PARTICLE particles[PARTICLE_NUMBER];
/************************************************************************************************************************/
/**** 如果采用这个onMouse()函数的话，则可以画出鼠标拖动矩形框的4种情形 ****/
/************************************************************************************************************************/
void onMouse(int event,int x,int y,int,void*)
{
//Point origin;//不能在这个地方进行定义，因为这是基于消息响应的函数，执行完后origin就释放了，所以达不到效果。
if(select_flag)
{
select.x=MIN(origin.x,x);//不一定要等鼠标弹起才计算矩形框，而应该在鼠标按下开始到弹起这段时间实时计算所选矩形框
select.y=MIN(origin.y,y);
select.width=abs(x-origin.x);//算矩形宽度和高度
select.height=abs(y-origin.y);
select&=Rect(0,0,frame.cols,frame.rows);//保证所选矩形框在视频显示区域之内
// rectangle(frame,select,Scalar(0,0,255),3,8,0);//显示手动选择的矩形框
}
if(event==CV_EVENT_LBUTTONDOWN)
{
select_flag=true;//鼠标按下的标志赋真值
tracking=false;
select_show=true;
after_select_frames=0;//还没开始选择，或者重新开始选择，计数为0
origin=Point(x,y);//保存下来单击是捕捉到的点
select=Rect(x,y,0,0);//这里一定要初始化，因为在opencv中Rect矩形框类内的点是包含左上角那个点的，但是不含右下角那个点。
}
else if(event==CV_EVENT_LBUTTONUP)
{
select_flag=false;
tracking=true;
select_show=false;
after_select_frames=1;//选择完后的那一帧当做第1帧
}
}
/****粒子权值降序排列函数****/
int particle_decrease(const void *p1,const void *p2)
{
PARTICLE* _p1=(PARTICLE*)p1;
PARTICLE* _p2=(PARTICLE*)p2;
if(_p1->weight<_p2->weight)
return 1;
else if(_p1->weight>_p2->weight)
return -1;
return 0;//相等的情况下返回0
}
int main(int argc, unsigned char* argv[])
{
char c;
Mat target_img,track_img;
Mat target_hist,track_hist;
PARTICLE *pParticle;
/***打开摄像头****/
VideoCapture cam(0);
if (!cam.isOpened())
return -1;
/****读取一帧图像****/
cam>>frame;
if(frame.empty())
return -1;
VideoWriter output_dst( "demo.avi", CV_FOURCC('M', 'J', 'P', 'G'), 10, frame.size(), 1 );
/****建立窗口****/
namedWindow("camera",1);//显示视频原图像的窗口
/****捕捉鼠标****/
setMouseCallback("camera",onMouse,0);
while(1)
{
/****读取一帧图像****/
cam>>frame;
if(frame.empty())
return -1;
/****将rgb空间转换为hsv空间****/
cvtColor(frame,hsv,CV_BGR2HSV);
if(tracking)
{
if(1==after_select_frames)//选择完目标区域后
{
/****计算目标模板的直方图特征****/
target_img=Mat(hsv,select);//在此之前先定义好target_img,然后这样赋值也行，要学会Mat的这个操作
calcHist(&target_img,1,channels,Mat(),target_hist,3,hist_size,ranges);
normalize(target_hist,target_hist);
/****初始化目标粒子****/
pParticle=particles;//指针初始化指向particles数组
for(int x=0;x<PARTICLE_NUMBER;x++)
{
pParticle->x=cvRound(select.x+0.5*select.width);//选定目标矩形框中心为初始粒子窗口中心
pParticle->y=cvRound(select.y+0.5*select.height);
pParticle->orix=pParticle->x;//粒子的原始坐标为选定矩形框(即目标)的中心
pParticle->oriy=pParticle->y;
pParticle->prex=pParticle->x;//更新上一次的粒子位置
pParticle->prey=pParticle->y;
pParticle->rect=select;
pParticle->prescale=1;
pParticle->scale=1;
pParticle->hist=target_hist;
pParticle->weight=0;
pParticle++;
}
}
else if(2==after_select_frames)//从第二帧开始就可以开始跟踪了
{
double sum=0.0;
pParticle=particles;
RNG rng;//随机数产生器
/****更新粒子结构体的大部分参数****/
for(int i=0;i<PARTICLE_NUMBER;i++)
{
int x,y;
int xpre,ypre;
double s,pres;
xpre=pParticle->x;
ypre=pParticle->y;
pres=pParticle->scale;
/****更新粒子的矩形区域即粒子中心****/
x=cvRound(A1*(pParticle->x-pParticle->orix)+A2*(pParticle->prex-pParticle->orix)+
B0*rng.gaussian(sigmax)+pParticle->orix);
pParticle->x=max(0,min(x,frame.cols-1));
y=cvRound(A1*(pParticle->y-pParticle->oriy)+A2*(pParticle->prey-pParticle->oriy)+
B0*rng.gaussian(sigmay)+pParticle->oriy);
pParticle->y=max(0,min(y,frame.rows-1));
s=A1*(pParticle->scale-1)+A2*(pParticle->prescale-1)+B0*(rng.gaussian(sigmas))+1.0;
pParticle->scale=max(1.0,min(s,3.0));
pParticle->prex=xpre;
pParticle->prey=ypre;
pParticle->prescale=pres;
// pParticle->orix=pParticle->orix;
// pParticle->oriy=pParticle->oriy;
//注意在c语言中，x-1.0，如果x是int型，则这句语法有错误,但如果前面加了cvRound(x-0.5)则是正确的
pParticle->rect.x=max(0,min(cvRound(pParticle->x-0.5*pParticle->scale*pParticle->rect.width),frame.cols));
pParticle->rect.y=max(0,min(cvRound(pParticle->y-0.5*pParticle->scale*pParticle->rect.height),frame.rows));
pParticle->rect.width=min(cvRound(pParticle->rect.width),frame.cols-pParticle->rect.x);
pParticle->rect.height=min(cvRound(pParticle->rect.height),frame.rows-pParticle->rect.y);
// pParticle->rect.width=min(cvRound(pParticle->scale*pParticle->rect.width),frame.cols-pParticle->rect.x);
// pParticle->rect.height=min(cvRound(pParticle->scale*pParticle->rect.height),frame.rows-pParticle->rect.y);
/****计算粒子区域的新的直方图特征****/
track_img=Mat(hsv,pParticle->rect);
calcHist(&track_img,1,channels,Mat(),track_hist,3,hist_size,ranges);
normalize(track_hist,track_hist);
/****更新粒子的权值****/
// pParticle->weight=compareHist(target_hist,track_hist,CV_COMP_INTERSECT);
//采用巴氏系数计算相似度,永远与最开始的那一目标帧相比较
pParticle->weight=1.0-compareHist(target_hist,track_hist,CV_COMP_BHATTACHARYYA);
/****累加粒子权值****/
sum+=pParticle->weight;
pParticle++;
}
/****归一化粒子权重****/
pParticle=particles;
for(int i=0;i<PARTICLE_NUMBER;i++)
{
pParticle->weight/=sum;
pParticle++;
}
/****根据粒子的权值降序排列****/
pParticle=particles;
qsort(pParticle,PARTICLE_NUMBER,sizeof(PARTICLE),&particle_decrease);
/****根据粒子权重重采样粒子****/
PARTICLE newParticle[PARTICLE_NUMBER];
int np=0,k=0;
for(int i=0;i<PARTICLE_NUMBER;i++)
{
np=cvRound(pParticle->weight*PARTICLE_NUMBER);
for(int j=0;j<np;j++)
{
newParticle[k++]=particles[i];
if(k==PARTICLE_NUMBER)
goto EXITOUT;
}
}
while(k<PARTICLE_NUMBER)
newParticle[k++]=particles[0];
EXITOUT:
for(int i=0;i<PARTICLE_NUMBER;i++)
particles[i]=newParticle[i];
}//end else
//????????这个排序很慢，粒子数一多就卡
// qsort(pParticle,PARTICLE_NUMBER,sizeof(PARTICLE),&particle_decrease);
/****计算粒子期望，采用所有粒子位置的期望值做为跟踪结果****/
/*Rect_<double> rectTrackingTemp(0.0,0.0,0.0,0.0);
pParticle=particles;
for(int i=0;i<PARTICLE_NUMBER;i++)
{
rectTrackingTemp.x+=pParticle->rect.x*pParticle->weight;
rectTrackingTemp.y+=pParticle->rect.y*pParticle->weight;
rectTrackingTemp.width+=pParticle->rect.width*pParticle->weight;
rectTrackingTemp.height+=pParticle->rect.height*pParticle->weight;
pParticle++;
}*/
/****计算最大权重目标的期望位置，作为跟踪结果****/
Rect rectTrackingTemp(0,0,0,0);
pParticle=particles;
rectTrackingTemp.x=pParticle->x-0.5*pParticle->rect.width;
rectTrackingTemp.y=pParticle->y-0.5*pParticle->rect.height;
rectTrackingTemp.width=pParticle->rect.width;
rectTrackingTemp.height=pParticle->rect.height;
/****计算最大权重目标的期望位置，采用权值最大的1/4个粒子数作为跟踪结果****/
/*Rect rectTrackingTemp(0,0,0,0);
double weight_temp=0.0;
pParticle=particles;
for(int i=0;i<PARTICLE_NUMBER/4;i++)
{
weight_temp+=pParticle->weight;
pParticle++;
}
pParticle=particles;
for(int i=0;i<PARTICLE_NUMBER/4;i++)
{
pParticle->weight/=weight_temp;
pParticle++;
}
pParticle=particles;
for(int i=0;i<PARTICLE_NUMBER/4;i++)
{
rectTrackingTemp.x+=pParticle->rect.x*pParticle->weight;
rectTrackingTemp.y+=pParticle->rect.y*pParticle->weight;
rectTrackingTemp.width+=pParticle->rect.width*pParticle->weight;
rectTrackingTemp.height+=pParticle->rect.height*pParticle->weight;
pParticle++;
}*/
/****计算最大权重目标的期望位置，采用所有粒子数作为跟踪结果****/
/*Rect rectTrackingTemp(0,0,0,0);
pParticle=particles;
for(int i=0;i<PARTICLE_NUMBER;i++)
{
rectTrackingTemp.x+=cvRound(pParticle->rect.x*pParticle->weight);
rectTrackingTemp.y+=cvRound(pParticle->rect.y*pParticle->weight);
pParticle++;
}
pParticle=particles;
rectTrackingTemp.width = pParticle->rect.width;
rectTrackingTemp.height = pParticle->rect.height;*/
//创建目标矩形区域
Rect tracking_rect(rectTrackingTemp);
pParticle=particles;
/****显示各粒子运动结果****/
for(int m=0;m<PARTICLE_NUMBER;m++)
{
rectangle(frame,pParticle->rect,Scalar(255,0,0),1,8,0);
pParticle++;
}
/****显示跟踪结果****/
rectangle(frame,tracking_rect,Scalar(0,0,255),3,8,0);
after_select_frames++;//总循环每循环一次，计数加1
if(after_select_frames>2)//防止跟踪太长，after_select_frames计数溢出
after_select_frames=2;
}
if(select_show)
rectangle(frame,select,Scalar(0,0,255),3,8,0);//显示手动选择的矩形框
output_dst<<frame;
//显示视频图片到窗口
imshow("camera",frame);
// select.zeros();
//键盘响应
c=(char)waitKey(20);
if(27==c)//ESC键
return -1;
}
return 0;
}结果：
扫一扫
4248
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
