基于粒子滤波的目标跟踪——解读Rob Hess的好文章

http://www.cnblogs.com/milier-otw/archive/2013/03/07.html

　　“一直都觉得粒子滤波是个挺牛的东西，每次试图看文献都被复杂的数学符号搞得看不下去。一个偶然的机会发现了Rob Hess(http://web.engr.oregonstate.edu/~hess/)实现的这个粒子滤波。从代码入手，一下子就明白了粒子滤波的原理。根据维基百科上对粒子滤波的介绍(http://en.wikipedia.org/wiki/Particle_filter)，粒子滤波其实有很多变种，Rob Hess实现的这种应该是最基本的一种，Sampling Importance Resampling (SIR)，根据重要性重采样。下面是我对粒子滤波实现物体跟踪的算法原理的粗浅理解

……”

详见如下网址：

http://www.cnblogs.com/yangyangcv/archive/2010/05/23/1742263.html

posted @ 2013-03-07 17:50 milier_otw 阅读(14) 评论(0)   编辑

GSL库在VC6.0上的配置

(转至http://blog.csdn.net/leowangzi/article/details/5357163)

这是我在VC6.0上配置Gsl的过程，希望对学习gsl的朋友有所帮助。

一、GSL介绍
GNU科学计算函数库GSL（GNU Scientific Library）是一个强大的C/C++数值计算函数库，它是一个自由软件，是GNU项目软件的一个部分，遵循GPL协议。GSL是一个为C和C++程序员提供的科学数值运算库。该科学计算库异常强大，函数库提供了大量的数值计算程序，如随机函数、特殊函数和拟合函数等等，整个函数库大约有1000多个函数，几乎涵盖了科学计算的各个方面。提供了如下方面的支持：
Complex Numbers Roots of Polynomials Special Functions
Vectors and Matrices Permutations Sorting
BLAS Support Linear Algebra Eigensystems
Fast Fourier Transforms Quadrature Random Numbers
Quasi-Random Sequences Random Distributions Statistics
Histograms N-Tuples Monte Carlo Integration
Simulated Annealing Differential Equations Interpolation
Numerical Differentiation Chebyshev Approximation Series Acceleration
Discrete Hankel Transforms Root-Finding Minimization
Least-Squares Fitting Physical Constants IEEE Floating-Point
Discrete Wavelet Transforms Basis splines

该函数库的主页是：http://www.gnu.org/software/gsl/gsl.html。

不过遗憾的是原始GSL并不支持不支持windows平台，可所幸的是有人做了GSL在windows上的移植工作，详见http://gnuwin32.sourceforge.net/packages/gsl.htm，目前版本是1.8。
二、下载gsl
1、从http://gnuwin32.sourceforge.net/packages/gsl.htm下载Complete package, except sources和Sources两个exe文件。
2、从http://www6.in.tum.de/~kiss/WinGsl.htm下载WinGsl-Lib-1.4.02.zip。
三、安装
1、 首先安装从http://gnuwin32.sourceforge.net/packages/gsl.htm下载的两个文件gsl-1.8.exe和gsl-1.8-src.exe。
2、 解压WinGsl-Lib-1.4.02.zip到D盘下，即D:/WinGsl。
四、设置Visual C++ 6.0编译环境
1、将D:/WinGsl/bin中的WinGsl.dll和WinGslD.dll复制到D:/VC6.0/Bin；将整个Gsl目录复制到D:/VC6.0/Bin下；lib目录下的所有.lib文件全部复制到D:/VC6.0/Lib下。
2、新建一个工程用于测试，我新建的是Test_gsl。然后在该项目的project－settings－link，在object/library modules中加入你用到的库文件，如WinGsl.lib等，加入多个可以用空格隔开。
3、 在tools－options－directories中，将D:/WinGsl下的lib，gsl分别加入到库文件和头文件的搜索路径中。
五、测试Gsl函数库
在工程Test_gsl中添加源文件test_gsl，其代码如下：
#include
#include
int main()
{
std::cout << gsl_sf_gamma_inc(1.5,0.5) << std::endl;
std::cout << gsl_sf_gamma_inc_Q( 1.5, 0.5 ) << std::endl;
std::cout << gsl_sf_gamma_inc_P( 1.5, 0.5 ) << std::endl;
return NULL;
}
结果运行的时候出现以下错误：
　　LIBCD.lib(wincrt0.obj) : error LNK2001: unresolved external symbol _WinMain@16
　　Debug/Test_gsl.exe : fatal error LNK1120: 1 unresolved externals
然后我们将在project->settings->link->project options里，把project设置中的/subsysetm:Windows改为/subsystem:console，便可以编译并运行Test_gsl工程，说明Gsl函数库已经可以使用了。

posted @ 2013-03-07 17:43 milier_otw 阅读(14) 评论(0)   编辑

学习OpenMP

最近在看多核编程。简单来说，由于现在电脑CPU一般都有两个核，4核与8核的CPU也逐渐走入了寻常百姓家，传统的单线程编程方式难以发挥多核CPU的强大功能，于是多核编程应运而生。按照我的理解，多核编程可以认为是对多线程编程做了一定程度的抽象，提供一些简单的API，使得用户不必花费太多精力来了解多线程的底层知识，从而提高编程效率。这两天关注的多核编程的工具包括openMP和TBB。按照目前网上的讨论，TBB风头要盖过openMP，比如openCV过去是使用openMP的，但从2.3版本开始抛弃openMP，转向TBB。但我试下来，TBB还是比较复杂的，相比之下，openMP则非常容易上手。因为精力和时间有限，没办法花费太多时间去学习TBB，就在这里分享下这两天学到的openMP的一点知识，和大家共同讨论。

openMP支持的编程语言包括C语言、C++和Fortran，支持OpenMP的编译器包括Sun Studio，Intel Compiler，Microsoft Visual Studio，GCC。我使用的是Microsoft Visual Studio 2008，CPU为Intel i5 四核，首先讲一下在Microsoft Visual Studio 2008上openMP的配置。非常简单，总共分2步：

(1) 新建一个工程。这个不再多讲。

(2) 建立工程后，点击 菜单栏->Project->Properties，弹出菜单里，点击 Configuration Properties->C/C++->Language->OpenMP Support，在下拉菜单里选择Yes。

至此配置结束。下面我们通过一个小例子来说明openMP的易用性。这个例子是 有一个简单的test()函数，然后在main()里，用一个for循环把这个test()函数跑8遍。

 1 #include <iostream>
 2 #include <time.h>
 3 void test()
 4 {
 5     int a = 0;
 6     for (int i=0;i<100000000;i++)
 7         a++;
 8 }
 9 int main()
10 {
11     clock_t t1 = clock();
12     for (int i=0;i<8;i++)
13         test();
14     clock_t t2 = clock();
15     std::cout<<"time: "<<t2-t1<<std::endl;
16 }

编译运行后，打印出来的耗时为：1.971秒。下面我们用一句话把上面代码变成多核运行。

 1 #include <iostream>
 2 #include <time.h>
 3 void test()
 4 {
 5     int a = 0;
 6     for (int i=0;i<100000000;i++)
 7         a++;
 8 }
 9 int main()
10 {
11     clock_t t1 = clock();
12     #pragma omp parallel for
13     for (int i=0;i<8;i++)
14         test();
15     clock_t t2 = clock();
16     std::cout<<"time: "<<t2-t1<<std::endl;
17 }

编译运行后，打印出来的耗时为：0.546秒，几乎为上面时间的1/4。

由此我们可以看到openMP的简单易用。在上面的代码里，我们一没有额外include头文件，二没有额外link库文件，只是在for循环前加了一句#pragma omp parallel for。而且这段代码在单核机器上，或者编译器没有将openMP设为Yes的机器上编译也不会报错，将自动忽略#pragma这行代码，然后按照传统单核串行的方式编译运行！我们唯一要多做的一步，是从C:\Program Files\Microsoft Visual Studio 9.0\VC\redist\x86\Microsoft.VC90.OPENMP和C:\Program Files\Microsoft Visual Studio 9.0\VC\redist\Debug_NonRedist\x86\Microsoft.VC90.DebugOpenMP目录下分别拷贝vcomp90d.dll和vcomp90.dll文件到工程文件当前目录下。

对上面代码按照我的理解做个简单的剖析。

当编译器发现#pragma omp parallel for后，自动将下面的for循环分成N份，(N为电脑CPU核数)，然后把每份指派给一个核去执行，而且多核之间为并行执行。下面的代码验证了这种分析。

1 #include <iostream>
2 int main()
3 {
4 #pragma omp parallel for
5     for (int i=0;i<10;i++)
6         std::cout<<i<<std::endl;
7     return 0;
8 }

会发现控制台打印出了0 3 4 5 8 9 6 7 1 2。注意：因为每个核之间是并行执行，所以每次执行时打印出的顺序可能都是不一样的。

下面我们来了谈谈竞态条件(race condition)的问题，这是所有多线程编程最棘手的问题。该问题可表述为，当多个线程并行执行时，有可能多个线程同时对某变量进行了读写操作，从而导致不可预知的结果。比如下面的例子，对于包含10个整型元素的数组a，我们用for循环求它各元素之和，并将结果保存在变量sum里。

 1 #include <iostream>
 2 int main()
 3 {
 4     int sum = 0;
 5     int a[10] = {1,2,3,4,5,6,7,8,9,10};
 6 #pragma omp parallel for
 7     for (int i=0;i<10;i++)
 8         sum = sum + a[i];
 9     std::cout<<"sum: "<<sum<<std::endl;
10     return 0;
11 }

如果我们注释掉#pragma omp parallel for，让程序先按照传统串行的方式执行，很明显，sum = 55。但按照并行方式执行后，sum则会变成其他值，比如在某次运行过程中，sum = 49。其原因是，当某线程A执行sum = sum + a[i]的同时，另一线程B正好在更新sum，而此时A还在用旧的sum做累加，于是出现了错误。

那么用openMP怎么实现并行数组求和呢？下面我们先给出一个基本的解决方案。该方案的思想是，首先生成一个数组sumArray，其长度为并行执行的线程的个数(默认情况下，该个数等于CPU的核数)，在for循环里，让各个线程更新自己线程对应的sumArray里的元素，最后再将sumArray里的元素累加到sum里，代码如下

 1 #include <iostream>
 2 #include <omp.h>
 3 int main(){
 4     int sum = 0;
 5     int a[10] = {1,2,3,4,5,6,7,8,9,10};
 6     int coreNum = omp_get_num_procs();//获得处理器个数
 7     int* sumArray = new int[coreNum];//对应处理器个数，先生成一个数组
 8     for (int i=0;i<coreNum;i++)//将数组各元素初始化为0
 9         sumArray[i] = 0;
10 #pragma omp parallel for
11     for (int i=0;i<10;i++)
12     {
13         int k = omp_get_thread_num();//获得每个线程的ID
14         sumArray[k] = sumArray[k]+a[i];
15     }
16     for (int i = 0;i<coreNum;i++)
17         sum = sum + sumArray[i];
18     std::cout<<"sum: "<<sum<<std::endl;
19     return 0;
20 }

需要注意的是，在上面代码里，我们用omp_get_num_procs()函数来获取处理器个数，用omp_get_thread_num()函数来获得每个线程的ID，为了使用这两个函数，我们需要include <omp.h>。

上面的代码虽然达到了目的，但它产生了较多的额外操作，比如要先生成数组sumArray，最后还要用一个for循环将它的各元素累加起来，有没有更简便的方式呢？答案是有，openMP为我们提供了另一个工具，归约(reduction)，见下面代码：

 1 #include <iostream>
 2 int main(){
 3     int sum = 0;
 4     int a[10] = {1,2,3,4,5,6,7,8,9,10};
 5 #pragma omp parallel for reduction(+:sum)
 6     for (int i=0;i<10;i++)
 7         sum = sum + a[i];
 8     std::cout<<"sum: "<<sum<<std::endl;
 9     return 0;
10 }

上面代码里，我们在#pragma omp parallel for 后面加上了 reduction(+:sum)，它的意思是告诉编译器：下面的for循环你要分成多个线程跑，但每个线程都要保存变量sum的拷贝，循环结束后，所有线程把自己的sum累加起来作为最后的输出。

reduction虽然很方便，但它只支持一些基本操作，比如+,-,*,&,|,&&,||等。有些情况下，我们既要避免race condition，但涉及到的操作又超出了reduction的能力范围，应该怎么办呢？这就要用到openMP的另一个工具，critical。来看下面的例子，该例中我们求数组a的最大值，将结果保存在max里。

 1 #include <iostream>
 2 int main(){
 3     int max = 0;
 4     int a[10] = {11,2,33,49,113,20,321,250,689,16};
 5 #pragma omp parallel for
 6     for (int i=0;i<10;i++)
 7     {
 8         int temp = a[i];
9  #pragma omp critical 
10         {
11             if (temp > max)
12                 max = temp;
13         }
14     }
15     std::cout<<"max: "<<max<<std::endl;
16     return 0;
17 }

上例中，for循环还是被自动分成N份来并行执行，但我们用#pragma omp critical将 if (temp > max) max = temp 括了起来，它的意思是：各个线程还是并行执行for里面的语句，但当你们执行到critical里面时，要注意有没有其他线程正在里面执行，如果有的话，要等其他线程执行完再进去执行。这样就避免了race condition问题，但显而易见，它的执行速度会变低，因为可能存在线程等待的情况。
有了以上基本知识，对我来说做很多事情都足够了。下面我们来看一个具体的应用例，从硬盘读入两幅图像，对这两幅图像分别提取特征点，特征点匹配，最后将图像与匹配特征点画出来。理解该例子需要一些图像处理的基本知识，我不在此详细介绍。另外，编译该例需要opencv，我用的版本是2.3.1，关于opencv的安装与配置也不在此介绍。我们首先来看传统串行编程的方式。

 1 #include "opencv2/highgui/highgui.hpp"
 2 #include "opencv2/features2d/features2d.hpp"
 3 #include <iostream>
 4 #include <omp.h>
 5 int main( ){
 6     cv::SurfFeatureDetector detector( 400 );    
 7     cv::SurfDescriptorExtractor extractor;
 8     cv::BruteForceMatcher<cv::L2<float> > matcher;
 9     std::vector< cv::DMatch > matches;
10     cv::Mat im0,im1;
11     std::vector<cv::KeyPoint> keypoints0,keypoints1;
12     cv::Mat descriptors0, descriptors1;
13     double t1 = omp_get_wtime( );
14     //先处理第一幅图像
15     im0 = cv::imread("rgb0.jpg", CV_LOAD_IMAGE_GRAYSCALE );
16     detector.detect( im0, keypoints0);
17     extractor.compute( im0,keypoints0,descriptors0);
18     std::cout<<"find "<<keypoints0.size()<<"keypoints in im0"<<std::endl;
19     //再处理第二幅图像
20     im1 = cv::imread("rgb1.jpg", CV_LOAD_IMAGE_GRAYSCALE );
21     detector.detect( im1, keypoints1);
22     extractor.compute( im1,keypoints1,descriptors1);
23     std::cout<<"find "<<keypoints1.size()<<"keypoints in im1"<<std::endl;
24     double t2 = omp_get_wtime( );
25     std::cout<<"time: "<<t2-t1<<std::endl;
26     matcher.match( descriptors0, descriptors1, matches );
27     cv::Mat img_matches;
28     cv::drawMatches( im0, keypoints0, im1, keypoints1, matches, img_matches ); 
29     cv::namedWindow("Matches",CV_WINDOW_AUTOSIZE);
30     cv::imshow( "Matches", img_matches );
31     cv::waitKey(0);
32     return 1;
33 }

很明显，读入图像，提取特征点与特征描述子这部分可以改为并行执行，修改如下：

 1 #include "opencv2/highgui/highgui.hpp"
 2 #include "opencv2/features2d/features2d.hpp"
 3 #include <iostream>
 4 #include <vector>
 5 #include <omp.h>
 6 int main( ){
 7     int imNum = 2;
 8     std::vector<cv::Mat> imVec(imNum);
 9     std::vector<std::vector<cv::KeyPoint>>keypointVec(imNum);
10     std::vector<cv::Mat> descriptorsVec(imNum);
11     cv::SurfFeatureDetector detector( 400 );    cv::SurfDescriptorExtractor extractor;
12     cv::BruteForceMatcher<cv::L2<float> > matcher;
13     std::vector< cv::DMatch > matches;
14     char filename[100];
15     double t1 = omp_get_wtime( );
16 #pragma omp parallel for
17     for (int i=0;i<imNum;i++){
18         sprintf(filename,"rgb%d.jpg",i);
19         imVec[i] = cv::imread( filename, CV_LOAD_IMAGE_GRAYSCALE );
20         detector.detect( imVec[i], keypointVec[i] );
21         extractor.compute( imVec[i],keypointVec[i],descriptorsVec[i]);
22         std::cout<<"find "<<keypointVec[i].size()<<"keypoints in im"<<i<<std::endl;
23     }
24     double t2 = omp_get_wtime( );
25     std::cout<<"time: "<<t2-t1<<std::endl;
26     matcher.match( descriptorsVec[0], descriptorsVec[1], matches );
27     cv::Mat img_matches;
28     cv::drawMatches( imVec[0], keypointVec[0], imVec[1], keypointVec[1], matches, img_matches ); 
29     cv::namedWindow("Matches",CV_WINDOW_AUTOSIZE);
30     cv::imshow( "Matches", img_matches );
31     cv::waitKey(0);
32     return 1;
33 }

两种执行方式做比较，时间为：2.343秒v.s. 1.2441秒

在上面代码中，为了改成适合#pragma omp parallel for执行的方式，我们用了STL的vector来分别存放两幅图像、特征点与特征描述子，但在某些情况下，变量可能不适合放在vector里，此时应该怎么办呢？这就要用到openMP的另一个工具，section，代码如下：

 1 #include "opencv2/highgui/highgui.hpp"
 2 #include "opencv2/features2d/features2d.hpp"
 3 #include <iostream>
 4 #include <omp.h>
 5 int main( ){
 6     cv::SurfFeatureDetector detector( 400 );    cv::SurfDescriptorExtractor extractor;
 7     cv::BruteForceMatcher<cv::L2<float> > matcher;
 8     std::vector< cv::DMatch > matches;
 9     cv::Mat im0,im1;
10     std::vector<cv::KeyPoint> keypoints0,keypoints1;
11     cv::Mat descriptors0, descriptors1;
12     double t1 = omp_get_wtime( );
13 #pragma omp parallel sections
14     {
15 #pragma omp section
16         {
17             std::cout<<"processing im0"<<std::endl;
18             im0 = cv::imread("rgb0.jpg", CV_LOAD_IMAGE_GRAYSCALE );
19             detector.detect( im0, keypoints0);
20             extractor.compute( im0,keypoints0,descriptors0);
21             std::cout<<"find "<<keypoints0.size()<<"keypoints in im0"<<std::endl;
22         }
23 #pragma omp section
24         {
25             std::cout<<"processing im1"<<std::endl;
26             im1 = cv::imread("rgb1.jpg", CV_LOAD_IMAGE_GRAYSCALE );
27             detector.detect( im1, keypoints1);
28             extractor.compute( im1,keypoints1,descriptors1);
29             std::cout<<"find "<<keypoints1.size()<<"keypoints in im1"<<std::endl;
30         }
31     }
32     double t2 = omp_get_wtime( );
33     std::cout<<"time: "<<t2-t1<<std::endl;
34     matcher.match( descriptors0, descriptors1, matches );
35     cv::Mat img_matches;
36     cv::drawMatches( im0, keypoints0, im1, keypoints1, matches, img_matches ); 
37     cv::namedWindow("Matches",CV_WINDOW_AUTOSIZE);
38     cv::imshow( "Matches", img_matches );
39     cv::waitKey(0);
40     return 1;
41 }

上面代码中，我们首先用#pragma omp parallel sections将要并行执行的内容括起来，在它里面，用了两个#pragma omp section，每个里面执行了图像读取、特征点与特征描述子提取。将其简化为伪代码形式即为：

 1 #pragma omp parallel sections
 2 {
 3     #pragma omp section
 4     {
 5         function1();
 6     }
 7 　　#pragma omp section
 8     {
 9         function2();
10     }
11 }

意思是：parallel sections里面的内容要并行执行，具体分工上，每个线程执行其中的一个section，如果section数大于线程数，那么就等某线程执行完它的section后，再继续执行剩下的section。在时间上，这种方式与人为用vector构造for循环的方式差不多，但无疑该种方式更方便，而且在单核机器上或没有开启openMP的编译器上，该种方式不需任何改动即可正确编译，并按照单核串行方式执行。

以上分享了这两天关于openMP的一点学习体会，其中难免有错误，欢迎指正。另外的一点疑问是，看到各种openMP教程里经常用到private,shared等来修饰变量，这些修饰符的意义和作用我大致明白，但在我上面所有例子中，不加这些修饰符似乎并不影响运行结果，不知道这里面有哪些讲究。