基于粒子滤波的物体跟踪
2014-09-20 15:10
204 查看
一直都觉得粒子滤波是个挺牛的东西，每次试图看文献都被复杂的数学符号搞得看不下去。一个偶然的机会发现了Rob
Hess(http://web.engr.oregonstate.edu/~hess/)实现的这个粒子滤波。从代码入手，一下子就明白了粒子滤波的原理。根据维基百科上对粒子滤波的介绍(http://en.wikipedia.org/wiki/Particle_filter)，粒子滤波其实有很多变种，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的思想也是(比如用在最简单的直线拟合上)，既然我不知道直线方程是什么，那我就随机的取两个点先算个直线出来，然后再看有多少点符合我的这条直线。哪条直线能获得最多的点的支持，哪条直线就是目标直线。想法非常简单，但效果很好。
代码和解释如下：
Hess(http://web.engr.oregonstate.edu/~hess/)实现的这个粒子滤波。从代码入手，一下子就明白了粒子滤波的原理。根据维基百科上对粒子滤波的介绍(http://en.wikipedia.org/wiki/Particle_filter)，粒子滤波其实有很多变种，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的思想也是(比如用在最简单的直线拟合上)，既然我不知道直线方程是什么，那我就随机的取两个点先算个直线出来，然后再看有多少点符合我的这条直线。哪条直线能获得最多的点的支持，哪条直线就是目标直线。想法非常简单，但效果很好。
代码和解释如下：
1 // particle_demo.cpp : 定义控制台应用程序的入口点。
2 //
3
4 #include "stdafx.h"
5 /************************************************************************/
6 /*
7 Description: 基本的粒子滤波目标跟踪
8 Author: Yang Xian
9 Email: yang_xian521@163.com
10 Version: 2011-11-2
11 History:
12 */
13 /************************************************************************/
14 #include <iostream> // for standard I/O
15 #include <string> // for strings
16 #include <iomanip> // for controlling float print precision
17 #include <sstream> // string to number conversion
18
19 #include <opencv2/imgproc/imgproc.hpp>
20 #include <opencv2/core/core.hpp> // Basic OpenCV structures (cv::Mat, Scalar)
21 #include <opencv2/highgui/highgui.hpp> // OpenCV window I/O
22
23 using namespace cv;
24 using namespace std;
25
26 // 以下这些参数对结果影响很大，而且也会根据视频内容，会对结果有很大的影响
27 const int PARTICLE_NUM = 25; // 粒子个数
28 // 粒子放入的相关区域
29 const double A1 = 2.0;
30 const double A2 = -1.0;
31 const double B0 = 1.0;
32 // 高斯随机数sigma参数
33 const double SIGMA_X = 1.0;
34 const double SIGMA_Y = 0.5;
35 const double SIGMA_SCALE = 0.001;
36
37 // 粒子结构体
38 typedef struct particle {
39 double x; // 当前x坐标
40 double y; // 当前y坐标
41 double scale; // 窗口比例系数
42 double xPre; // x坐标预测位置
43 double yPre; // y坐标预测位置
44 double scalePre; // 窗口预测比例系数
45 double xOri; // 原始x坐标
46 double yOri; // 原始y坐标
47 Rect rect; // 原始区域大小
48 MatND hist; // 粒子区域的特征直方图
49 double weight; // 该粒子的权重
50 } PARTICLE;
51
52 Mat hsv; // hsv色彩空间的输入图像
53 Mat roiImage; // 目标区域
54 MatND roiHist; // 目标区域直方图
55 Mat img; // 输出的目标图像
56 PARTICLE particles[PARTICLE_NUM]; // 粒子
57
58 int nFrameNum = 0;
59
60 bool bSelectObject = false; // 区域选择标志
61 bool bTracking = false; // 开始跟踪标志
62 Point origin; // 鼠标按下时的点位置
63 Rect selection;// 感兴趣的区域大小
64
65 // 直方图相关参数，特征的选取也会对结果影响巨大
66 // Quantize the hue to 30 levels
67 // and the saturation to 32 levels
68 // value to 10 levels
69 int hbins = 180, sbins = 256, vbin = 10;
70 int histSize[] = {hbins, sbins, vbin};
71 // hue varies from 0 to 179, see cvtColor
72 float hranges[] = { 0, 180 };
73 // saturation varies from 0 (black-gray-white) to 255 (pure spectrum color)
74 float sranges[] = { 0, 256 };
75 // value varies from 0 (black-gray-white) to 255 (pure spectrum color)
76 float vranges[] = { 0, 256 };
77 const float* ranges[] = {hranges, sranges, vranges};
78 // we compute the histogram from the 0-th and 1-st channels
79 int channels[] = {0, 1, 2};
80
81 // 鼠标响应函数，得到选择的区域，保存在selection
82 void onMouse(int event, int x, int y, int, void*)
83 {
84 if( bSelectObject )
85 {
86 selection.x = MIN(x, origin.x);
87 selection.y = MIN(y, origin.y);
88 selection.width = std::abs(x - origin.x);
89 selection.height = std::abs(y - origin.y);
90
91 selection &= Rect(0, 0, img.cols, img.rows);
92 }
93
94 switch (event)
95 {
96 case CV_EVENT_LBUTTONDOWN:
97 origin = Point(x,y);
98 selection = Rect(x,y,0,0);
99 bSelectObject = true;
100 bTracking = false;
101 break;
102 case CV_EVENT_LBUTTONUP:
103 bSelectObject = false;
104 // if( selection.width > 0 && selection.height > 0 )
105 bTracking = true;
106 nFrameNum = 0;
107 break;
108 }
109 }
110
111 // 快速排序算法排序函数
112 int particle_cmp(const void* p1,const void* p2)
113 {
114 PARTICLE* _p1 = (PARTICLE*)p1;
115 PARTICLE* _p2 = (PARTICLE*)p2;
116
117 if(_p1->weight < _p2->weight)
118 return 1; //按照权重降序排序
119 if(_p1->weight > _p2->weight)
120 return -1;
121 return 0;
122 }
123
124 const char* keys =
125 {
126 "{1| | 0 | camera number}"
127 };
128
129 int main(int argc, const char **argv)//这里char **argv前必须用const，why？
130 {
131 int delay = 30; // 控制播放速度
132 char c; // 键值
133
134 /*读取avi文件*/
135 VideoCapture captRefrnc("IndoorGTTest1.avi"); // 视频文件
136
137 /*打开摄像头*/
138 //VideoCapture captRefrnc;
139 //CommandLineParser parser(argc, argv, keys);//命令解析器函数
140 //int camNum = parser.get<int>("1");
141 //captRefrnc.open(camNum);
142
143 if ( !captRefrnc.isOpened())
144 {
145 return -1;
146 }
147
148 // Windows
149 const char* WIN_RESULT = "Result";
150 namedWindow(WIN_RESULT, CV_WINDOW_AUTOSIZE);
151 //namedWindow(WIN_RESULT, 0);
152 // namedWindow("Result",0);
153 // 鼠标响应函数
154 //setMouseCallback(WIN_RESULT, onMouse, 0);
155 setMouseCallback("Result", onMouse, 0);
156
157 Mat frame; //视频的每一帧图像
158
159 bool paused = false;
160 PARTICLE * pParticles = particles;//particles为可装PARTICLE_NUM个PARTICLE结构体的数组，所以pParticles为指向其数组的指针
161
162 while(true) //Show the image captured in the window and repeat
163 {
164 if(!paused)
165 {
166 captRefrnc >> frame;
167 if(frame.empty())
168 break;
169 }
170
171 frame.copyTo(img); // 接下来的操作都是对src的
172
173
174 // 选择目标后进行跟踪
175 if (bTracking == true)//鼠标操作选完后
176 {
177 if(!paused)
178 {
179 nFrameNum++;//帧数计数器
180 cvtColor(img, hsv, CV_BGR2HSV);
181 Mat roiImage(hsv, selection); // 目标区域，这个构造函数第二个参数表示截取selection部分
182
183 if (nFrameNum == 1) //选择目标后的第一帧需要初始化
184 {
185 // step 1: 提取目标区域特征，难道其目标特征就是其色调的直方图分布？
186 calcHist(&roiImage, 1, channels, Mat(), roiHist, 3, histSize, ranges);
187 normalize(roiHist, roiHist); // 归一化L2
188
189 // step 2: 初始化particle
190 pParticles = particles;
191 for (int i=0; i<PARTICLE_NUM; i++)
192 {
193 pParticles->x = selection.x + 0.5 * selection.width;
194 pParticles->y = selection.y + 0.5 * selection.height;
195 pParticles->xPre = pParticles->x;
196 pParticles->yPre = pParticles->y;
197 pParticles->xOri = pParticles->x;
198 pParticles->yOri = pParticles->y;
199 pParticles->rect = selection;
200 pParticles->scale = 1.0;
201 pParticles->scalePre = 1.0;
202 pParticles->hist = roiHist;
203 pParticles->weight = 0;
204 pParticles++;
205 }
206 }
207 else //不是第一帧
208 {
209 pParticles = particles;
210 RNG rng;//随机数序列产生器
211 for (int i=0; i<PARTICLE_NUM; i++)
212 {
213 // step 3: 求particle的transition，粒子结构中的参数全部更新过
214 double x, y, s;
215
216 pParticles->xPre = pParticles->x;
217 pParticles->yPre = pParticles->y;
218 pParticles->scalePre = pParticles->scale;
219
220 x = A1 * (pParticles->x - pParticles->xOri) + A2 * (pParticles->xPre - pParticles->xOri) +
221 B0 * rng.gaussian(SIGMA_X) + pParticles->xOri;//以当前点为中心产生高斯分布的粒子
222 pParticles->x = std::max(0.0, std::min(x, img.cols-1.0));//其实就是x，只是考虑了边界在内而已
223
224
225 y = A1 * (pParticles->y - pParticles->yOri) + A2 * (pParticles->yPre - pParticles->yOri) +
226 B0 * rng.gaussian(SIGMA_Y) + pParticles->yOri;
227 pParticles->y = std::max(0.0, std::min(y, img.rows-1.0));
228
229 s = A1 * (pParticles->scale - 1.0) + A2 * (pParticles->scalePre - 1.0) +
230 B0 * rng.gaussian(SIGMA_SCALE) + 1.0;
231 pParticles->scale = std::max(0.1, std::min(s, 3.0));
232 // rect参数有待考证
233 pParticles->rect.x = std::max(0, std::min(cvRound(pParticles->x - 0.5 * pParticles->rect.width * pParticles->scale), img.cols-1)); // 0 <= x <= img.rows-1
234 pParticles->rect.y = std::max(0, std::min(cvRound(pParticles->y - 0.5 * pParticles->rect.height * pParticles->scale), img.rows-1)); // 0 <= y <= img.cols-1
235 pParticles->rect.width = std::min(cvRound(pParticles->rect.width * pParticles->scale), img.cols - pParticles->rect.x);
236 pParticles->rect.height = std::min(cvRound(pParticles->rect.height * pParticles->scale), img.rows - pParticles->rect.y);
237 // Ori参数不改变
238
239 // step 4: 求particle区域的特征直方图
240 Mat imgParticle(img, pParticles->rect);
241 calcHist(&imgParticle, 1, channels, Mat(), pParticles->hist, 3, histSize, ranges);
242 normalize(pParticles->hist, pParticles->hist); // 归一化L2
243
244 // step 5: 特征的比对,更新particle权重
245 //compareHist()函数返回2个直方图之间的相似度，因为参数为CV_COMP_INTERSECT，所以返回的是最小直方图值之和
246 pParticles->weight = compareHist(roiHist, pParticles->hist, CV_COMP_INTERSECT);//其差值直接作为其权值
247
248 pParticles++;
249 }
250
251 // step 6: 归一化粒子权重
252 double sum = 0.0;
253 int i;
254
255 pParticles = particles;
256 for(i=0; i<PARTICLE_NUM; i++)
257 {
258 sum += pParticles->weight;
259 pParticles++;
260 }
261 pParticles = particles;
262 for(i=0; i<PARTICLE_NUM; i++)
263 {
264 pParticles->weight /= sum;
265 pParticles++;
266 }
267
268 // step 7: resample根据粒子的权重的后验概率分布重新采样
269 pParticles = particles;
270 // PARTICLE* newParticles = new PARTICLE[sizeof(PARTICLE) * PARTICLE_NUM];
271 PARTICLE newParticles[PARTICLE_NUM];
272 int np, k = 0;
273
274 //qsort()函数为对数组进行快速排序,其中第4个参数表示的是排序是升序还是降序
275 qsort(pParticles, PARTICLE_NUM, sizeof(PARTICLE), &particle_cmp);//这里采用的是降序排列
276 for(int i=0; i<PARTICLE_NUM; i++)
277 {
278 np = cvRound(particles[i].weight * PARTICLE_NUM);//权值高的优先重采样
279 for(int j=0; j<np; j++)
280 {
281 newParticles[k++] = particles[i];
282 if(k == PARTICLE_NUM)//重采样后达到了个数要求则直接跳出
283 goto EXITOUT;
284 }
285 }
286 while(k < PARTICLE_NUM)
287 {
288 newParticles[k++] = particles[0];//个数不够时，将权值最高的粒子重复给
289 }
290
291 EXITOUT:
292 for (int i=0; i<PARTICLE_NUM; i++)
293 {
294 particles[i] = newParticles[i];
295 }
296
297 }// end else
298
299 qsort(pParticles, PARTICLE_NUM, sizeof(PARTICLE), &particle_cmp);
300
301 // step 8: 计算粒子的期望，作为跟踪结果
302 Rect_<double> rectTrackingTemp(0.0, 0.0, 0.0, 0.0);
303 pParticles = particles;
304 for (int i=0; i<PARTICLE_NUM; i++)
305 {
306 rectTrackingTemp.x += pParticles->rect.x * pParticles->weight;//坐标加上权重的偏移值
307 rectTrackingTemp.y += pParticles->rect.y * pParticles->weight;
308 rectTrackingTemp.width += pParticles->rect.width * pParticles->weight;//宽度也要加上权值的偏移值
309 rectTrackingTemp.height += pParticles->rect.height * pParticles->weight;
310 pParticles++;
311 }
312 Rect rectTracking(rectTrackingTemp); // 跟踪结果
313
314 // 显示各粒子的运动
315 for (int i=0; i<PARTICLE_NUM; i++)
316 {
317 rectangle(img, particles[i].rect, Scalar(255,0,0));
318 }
319 // 显示跟踪结果
320 rectangle(img, rectTracking, Scalar(0,0,255), 3);
321
322 }
323 }// end Tracking
324
325 // imshow(WIN_SRC, frame);
326 imshow(WIN_RESULT, img);
327
328 c = (char)waitKey(delay);
329 if( c == 27 )
330 break;
331 switch(c)
332 {
333 case 'p'://暂停键
334 paused = !paused;
335 break;
336 default:
337 ;
338 }
339 }// end while
340 }
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- 基于粒子滤波的物体跟踪 Particle Filter Object Tracking
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪zz
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 转#基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪（应用代码，在LINUX下）
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 20170330_滤波_基于粒子滤波的物体跟踪
- [转]基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 基于粒子滤波的物体跟踪
- 【computer vision】基于粒子滤波的物体跟踪