在学习卡尔曼滤波器之前,首先看看为什么叫“卡尔曼”。跟其他著名的理论(例如傅立叶变换,泰勒级数等等)一样,卡尔曼也是一个人的名字,而跟他们不同的是,他是个现代人!
简单来说,卡尔曼滤波器是一个“optimalrecursivedataprocessingalgorithm(最优化自回归数据处理算法)”。对于解决很大部分的问题,他是最优,效率最高甚至是最有用的。他的广泛应用已经超过30年,包括机器人导航,控制,传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理,例如头脸识别,图像分割,图像边缘检测等等。
2.卡尔曼滤波器的介绍
(IntroductiontotheKalmanFilter)
为了可以更加容易的理解卡尔曼滤波器,这里会应用形象的描述方法来讲解,而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是,他的5条公式是其核心内容。结合现代的计算机,其实卡尔曼的程序相当的简单,只要你理解了他的那5条公式。
在介绍他的5条公式之前,先让我们来根据下面的例子一步一步的探索。
好了,现在对于某一分钟我们有两个有关于该房间的温度值:你根据经验的预测值(系统的预测值)和温度计的值(测量值)。下面我们要用这两个值结合他们各自的噪声来估算出房间的实际温度值。
假如我们要估算k时刻的是实际温度值。首先你要根据k-1时刻的温度值,来预测k时刻的温度。因为你相信温度是恒定的,所以你会得到k时刻的温度预测值是跟k-1时刻一样的,假设是23度,同时该值的高斯噪声的偏差是5度(5是这样得到的:如果k-1时刻估算出的最优温度值的偏差是3,你对自己预测的不确定度是4度,他们平方相加再开方,就是5)。然后,你从温度计那里得到了k时刻的温度值,假设是25度,同时该值的偏差是4度。
由于我们用于估算k时刻的实际温度有两个温度值,分别是23度和25度。究竟实际温度是多少呢?相信自己还是相信温度计呢?究竟相信谁多一点,我们可以用他们的covariance来判断。因为Kg^2=5^2/(5^2+4^2),所以Kg=0.78,我们可以估算出k时刻的实际温度值是:23+0.78*(25-23)=24.56度。可以看出,因为温度计的covariance比较小(比较相信温度计),所以估算出的最优温度值偏向温度计的值。
现在我们已经得到k时刻的最优温度值了,下一步就是要进入k+1时刻,进行新的最优估算。到现在为止,好像还没看到什么自回归的东西出现。对了,在进入k+1时刻之前,我们还要算出k时刻那个最优值(24.56度)的偏差。算法如下:((1-Kg)*5^2)^0.5=2.35。这里的5就是上面的k时刻你预测的那个23度温度值的偏差,得出的2.35就是进入k+1时刻以后k时刻估算出的最优温度值的偏差(对应于上面的3)。
就是这样,卡尔曼滤波器就不断的把covariance递归,从而估算出最优的温度值。他运行的很快,而且它只保留了上一时刻的covariance。上面的Kg,就是卡尔曼增益(KalmanGain)。他可以随不同的时刻而改变他自己的值,是不是很神奇!
下面就要言归正传,讨论真正工程系统上的卡尔曼。
3.卡尔曼滤波器算法
(TheKalmanFilterAlgorithm)
在这一部分,我们就来描述源于DrKalman的卡尔曼滤波器。下面的描述,会涉及一些基本的概念知识,包括概率(Probability),随即变量(RandomVariable),高斯或正态分配(GaussianDistribution)还有State-spaceModel等等。但对于卡尔曼滤波器的详细证明,这里不能一一描述。
首先,我们先要引入一个离散控制过程的系统。该系统可用一个线性随机微分方程(LinearStochasticDifferenceequation)来描述:
X(k)=AX(k-1)+BU(k)+W(k)
再加上系统的测量值:
Z(k)=HX(k)+V(k)
上两式子中,X(k)是k时刻的系统状态,U(k)是k时刻对系统的控制量。A和B是系统参数,对于多模型系统,他们为矩阵。Z(k)是k时刻的测量值,H是测量系统的参数,对于多测量系统,H为矩阵。W(k)和V(k)分别表示过程和测量的噪声。他们被假设成高斯白噪声(WhiteGaussianNoise),他们的covariance分别是Q,R(这里我们假设他们不随系统状态变化而变化)。
对于满足上面的条件(线性随机微分系统,过程和测量都是高斯白噪声),卡尔曼滤波器是最优的信息处理器。下面我们来用他们结合他们的covariances来估算系统的最优化输出(类似上一节那个温度的例子)。
首先我们要利用系统的过程模型,来预测下一状态的系统。假设现在的系统状态是k,根据系统的模型,可以基于系统的上一状态而预测出现在状态:
X(k|k-1)=AX(k-1|k-1)+BU(k)………..(1)
式(1)中,X(k|k-1)是利用上一状态预测的结果,X(k-1|k-1)是上一状态最优的结果,U(k)为现在状态的控制量,如果没有控制量,它可以为0。
到现在为止,我们的系统结果已经更新了,可是,对应于X(k|k-1)的covariance还没更新。我们用P表示covariance:
P(k|k-1)=AP(k-1|k-1)A’+Q………(2)
式(2)中,P(k|k-1)是X(k|k-1)对应的covariance,P(k-1|k-1)是X(k-1|k-1)对应的covariance,A’表示A的转置矩阵,Q是系统过程的covariance。式子1,2就是卡尔曼滤波器5个公式当中的前两个,也就是对系统的预测。
现在我们有了现在状态的预测结果,然后我们再收集现在状态的测量值。结合预测值和测量值,我们可以得到现在状态(k)的最优化估算值X(k|k):
X(k|k)=X(k|k-1)+Kg(k)(Z(k)-HX(k|k-1))………(3)
其中Kg为卡尔曼增益(KalmanGain):
Kg(k)=P(k|k-1)H’/(HP(k|k-1)H’+R)………(4)
到现在为止,我们已经得到了k状态下最优的估算值X(k|k)。但是为了要另卡尔曼滤波器不断的运行下去直到系统过程结束,我们还要更新k状态下X(k|k)的covariance:
P(k|k)=(I-Kg(k)H)P(k|k-1)………(5)
其中I为1的矩阵,对于单模型单测量,I=1。当系统进入k+1状态时,P(k|k)就是式子(2)的P(k-1|k-1)。这样,算法就可以自回归的运算下去。
卡尔曼滤波器的原理基本描述了,式子1,2,3,4和5就是他的5个基本公式。根据这5个公式,可以很容易的实现计算机的程序
CamShift算法简介CamShift算法,即"ContinuouslyApativeMean-Shift"算法,是一种运动跟踪算法。它主要通过视频图像中运动物体的颜色信息来达到跟踪的目的。我把这个算法分解成三个部分,便于理解:
BackProjection计算。
MeanShift算法
CamShift算法
1BackProjection计算计算BackProjection的步骤是这样的:
1.计算被跟踪目标的色彩直方图。在各种色彩空间中,只有HSI空间(或与HSI类似的色彩空间)中的H分量可以表示颜色信息。所以在具体的计算过程中,首先将其他的色彩空间的值转化到HSI空间,然后会其中的H分量做1D直方图计算。
2.根据获得的色彩直方图将原始图像转化成色彩概率分布图像,这个过程就被称作"BackProjection"。
在OpenCV中的直方图函数中,包含BackProjection的函数,函数原型是:
voidcvCalcBackProject(IplImage**img,CvArr**backproject,constCvHistogram*hist);
传递给这个函数的参数有三个:
1.IplImage**img:存放原始图像,输入。
2.CvArr**backproject:存放BackProjection结果,输出。
3.CvHistogram*hist:存放直方图,输入
下面就给出计算BackProjection的OpenCV代码。
1.准备一张只包含被跟踪目标的图片,将色彩空间转化到HSI空间,获得其中的H分量:
IplImage*target=cvLoadImage("target.bmp",-1);//装载图片
IplImage*target_hsv=cvCreateImage(cvGetSize(target),IPL_DEPTH_8U,3);
IplImage*target_hue=cvCreateImage(cvGetSize(target),IPL_DEPTH_8U,3);
cvCvtColor(target,target_hsv,CV_BGR2HSV);//转化到HSV空间
cvSplit(target_hsv,target_hue,NULL,NULL,NULL);//获得H分量
2.计算H分量的直方图,即1D直方图:
IplImage*h_plane=cvCreateImage(cvGetSize(target_hsv),IPL_DEPTH_8U,1);
inthist_size[]={255};//将H分量的值量化到[0,255]
float*ranges[]={{0,360}};//H分量的取值范围是[0,360)
CvHistogram*hist=cvCreateHist(1,hist_size,ranges,1);
cvCalcHist(&target_hue,hist,0,NULL);
在这里需要考虑H分量的取值范围的问题,H分量的取值范围是[0,360),这个取值范围的值不能用一个byte来表示,为了能用一个byte表示,需要将H值做适当的量化处理,在这里我们将H分量的范围量化到[0,255].
4.计算BackProjection:
IplImage*rawImage;
//----------------------------------------------
//getfromvideoframe,unsignedbyte,onechannel
IplImage*result=cvCreateImage(cvGetSize(rawImage),IPL_DEPTH_8U,1);
cvCalcBackProject(&rawImage,result,hist);
5.结果:result即为我们需要的.
2)MeanShift算法
这里来到了CamShift算法,OpenCV实现的第二部分,这一次重点讨论MeanShift算法。
在讨论MeanShift算法之前,首先讨论在2D概率分布图像中,如何计算某个区域的重心(MassCenter)的问题,重心可以通过以下公式来计算:
1.计算区域内0阶矩
for(inti=0;i for(intj=0;j M00+=I(i,j) 2.区域内1阶矩: { M10+=i*I(i,j); M01+=j*I(i,j); } 3.则MassCenter为: Xc=M10/M00;Yc=M01/M00 接下来,讨论MeanShift算法的具体步骤,MeanShift算法可以分为以下4步: 1.选择窗的大小和初始位置. 2.计算此时窗口内的MassCenter. 3.调整窗口的中心到MassCenter. 4.重复2和3,直到窗口中心"会聚",即每次窗口移动的距离小于一定的阈值。 在OpenCV中,提供MeanShift算法的函数,函数的原型是: intcvMeanShift(IplImage*imgprob,CvRectwindowIn, CvTermCriteriacriteria,CvConnectedComp*out); 需要的参数为: 1.IplImage*imgprob:2D概率分布图像,传入; 2.CvRectwindowIn:初始的窗口,传入; 3.CvTermCriteriacriteria:停止迭代的标准,传入; 4.CvConnectedComp*out:查询结果,传出。 (注:构造CvTermCriteria变量需要三个参数,一个是类型,另一个是迭代的最大次数,最后一个表示特定的阈值。例如可以这样构造criteria:criteria=cvTermCriteria(CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,10,0.1)。) 返回的参数: 1.int:迭代的次数。 实现代码:暂时缺 3)CamShift算法1.原理 在了解了MeanShift算法以后,我们将MeanShift算法扩展到连续图像序列(一般都是指视频图像序列),这样就形成了CamShift算法。CamShift算法的全称是"ContinuouslyApaptiveMean-SHIFT",它的基本思想是视频图像的所有帧作MeanShift运算,并将上一帧的结果(即SearchWindow的中心和大小)作为下一帧MeanShift算法的SearchWindow的初始值,如此迭代下去,就可以实现对目标的跟踪。整个算法的具体步骤分5步: Step1:将整个图像设为搜寻区域。 Step2:初始话SearchWindow的大小和位置。 Step3:计算SearchWindow内的彩色概率分布,此区域的大小比SearchWindow要稍微大一点。 Step4:运行MeanShift。获得SearchWindow新的位置和大小。 Step5:在下一帧视频图像中,用Step3获得的值初始化SearchWindow的位置和大小。跳转到Step3继续运行。 2.实现 在OpenCV中,有实现CamShift算法的函数,此函数的原型是: cvCamShift(IplImage*imgprob,CvRectwindowIn, CvTermCriteriacriteria, CvConnectedComp*out,CvBox2D*box=0); 其中: imgprob:色彩概率分布图像。 windowIn:SearchWindow的初始值。 Criteria:用来判断搜寻是否停止的一个标准。 out:保存运算结果,包括新的SearchWindow的位置和面积。 box:包含被跟踪物体的最小矩形。 说明: 1.在OpenCV4.0beta的目录中,有CamShift的例子。遗憾的是这个例子目标的跟踪是半自动的,即需要人手工选定一个目标。我正在努力尝试全自动的目标跟踪,希望可以和大家能在这方面与大家交流。 2006.03.15来自:CSDN选自Mailbomb的blog 程序员每天该做的事 1、总结自己一天任务的完成情况 最好的方式是写工作日志,把自己今天完成了什么事情,遇见了什么问题都记录下来,日后翻看好处多多 2、考虑自己明天应该做的主要工作 3、考虑自己一天工作中失误的地方,并想出避免下一次再犯的方法 出错不要紧,最重要的是不要重复犯相同的错误,那是愚蠢 4、考虑自己一天工作完成的质量和效率能否还能提高 一天只提高1%,365天你的效率就能提高多少倍你知道吗?(1+0.01)^365=37倍 5、看一个有用的新闻网站或读一张有用的报纸,了解业界动态 闭门造车是不行的,了解一下别人都在做什么,对自己能带来很多启示 6、记住一位同事的名字及其特点 你认识公司的所有同事吗?你了解他们吗? 7、清理自己的代码 今天完成的代码,把中间的调试信息,测试代码清理掉,按照编码风格整理好,注释都写好了吗? 8、清理自己的桌面 当日事当日毕,保持清洁干劲的桌面才能让你工作时不分心,程序员特别要把电脑的桌面清理干净 程序员每月该做的事 1、至少和一个同事一起吃饭或喝茶不光了解自己工作伙伴的工作,还要了解他们的生活 2、自我考核一次相对正式地考核自己一下,你对得起这个月的工资吗? 3、对你的同事考核一次你的同事表现怎么样?哪些人值得学习,哪些人需要帮助? 3、制定下月的计划,确定下月的工作重点 4、总结自己工作质量改进状况自己的质量提高了多少? 5、有针对性地对一项工作指标做深入地分析并得出改进的方案可以是对自己的,也可以是对公司的,一定要深入地分析后拿出自己的观点来。要想在老板面前说得上话,做的成事,工作上功夫要做足。 6、与老板沟通一次最好是面对面地沟通,好好表现一下自己,虚心听取老板的意见,更重要的是要了解老板当前关心的重点 程序员每年该做的事 1、年终总结每个公司都会做的事情,但你真正认真地总结过自己吗? 2、兑现给自己、给家人的承诺给老婆、儿子的新年礼物买了没有?给自己的呢? 3、下年度工作规划好好想想自己明年的发展目标,争取升职/加薪、跳槽还是自己出来干? 4、掌握一项新技术至少是一项,作为程序员一年要是一项新技术都学不到手,那就一定会被淘汰。掌握可不是看本书就行的,要真正懂得应用,最好你能够写一篇教程发表到你的blog 5、推出一种新产品可以是一个真正的产品,也可以只是一个类库,只要是你创造的东西就行,让别人使用它,也为世界作点贡献。当然如果真的很有价值,收点注册费也是应该的