2008年8月31日星期日

OpenCV 初级使用指导-1

矩阵操作

分配释放矩阵空间

  • 综述:
    • OpenCV有针对矩阵操作的C语言函数. 许多其他方法提供了更加方便的C++接口,其效率与OpenCV一样.
    • OpenCV将向量作为1维矩阵处理.
    • 矩阵按行存储,每行有4字节的校整.

  • 分配矩阵空间:
    CvMat* cvCreateMat(int rows, int cols, int type);
    

    type: 矩阵元素类型. 格式为CV_(S|U|F)C.
    例如: CV_8UC1 表示8位无符号单通道矩阵, CV_32SC2表示32位有符号双通道矩阵.

    例程:
    CvMat* M = cvCreateMat(4,4,CV_32FC1);

  • 释放矩阵空间:
    CvMat* M = cvCreateMat(4,4,CV_32FC1);
    
    cvReleaseMat(&M);

  • 复制矩阵:
    CvMat* M1 = cvCreateMat(4,4,CV_32FC1);
    
    CvMat* M2;
    M2=cvCloneMat(M1);

  • 初始化矩阵:
    double a[] = { 1,   2,   3,   4,
    
    5, 6, 7, 8,
    9, 10, 11, 12 };

    CvMat Ma=cvMat(3, 4, CV_64FC1, a);

    另一种方法:

    CvMat Ma;
    
    cvInitMatHeader(&Ma, 3, 4, CV_64FC1, a);

  • 初始化矩阵为单位阵:
    CvMat* M = cvCreateMat(4,4,CV_32FC1);
    
    cvSetIdentity(M); // 这里似乎有问题,不成功

存取矩阵元素

  • 假设需要存取一个2维浮点矩阵的第(i,j)个元素.

  • 间接存取矩阵元素:
    cvmSet(M,i,j,2.0); // Set M(i,j)
    
    t = cvmGet(M,i,j); // Get M(i,j)

  • 直接存取,假设使用4-字节校正:
    CvMat* M     = cvCreateMat(4,4,CV_32FC1);
    
    int n = M->cols;
    float *data = M->data.fl;

    data[i*n+j] = 3.0;

  • 直接存取,校正字节任意:
    CvMat* M     = cvCreateMat(4,4,CV_32FC1);
    
    int step = M->step/sizeof(float);
    float *data = M->data.fl;

    (data+i*step)[j] = 3.0;

  • 直接存取一个初始化的矩阵元素:
    double a[16];
    
    CvMat Ma = cvMat(3, 4, CV_64FC1, a);
    a[i*4+j] = 2.0; // Ma(i,j)=2.0;

矩阵/向量操作

  • 矩阵-矩阵操作:
    CvMat *Ma, *Mb, *Mc;
    
    cvAdd(Ma, Mb, Mc); // Ma+Mb -> Mc
    cvSub(Ma, Mb, Mc); // Ma-Mb -> Mc
    cvMatMul(Ma, Mb, Mc); // Ma*Mb -> Mc

  • 按元素的矩阵操作:
    CvMat *Ma, *Mb, *Mc;
    
    cvMul(Ma, Mb, Mc); // Ma.*Mb -> Mc
    cvDiv(Ma, Mb, Mc); // Ma./Mb -> Mc
    cvAddS(Ma, cvScalar(-10.0), Mc); // Ma.-10 -> Mc

  • 向量乘积:
    double va[] = {1, 2, 3};
    
    double vb[] = {0, 0, 1};
    double vc[3];

    CvMat Va=cvMat(3, 1, CV_64FC1, va);
    CvMat Vb=cvMat(3, 1, CV_64FC1, vb);
    CvMat Vc=cvMat(3, 1, CV_64FC1, vc);

    double res=cvDotProduct(&Va,&Vb); // 点乘: Va . Vb -> res
    cvCrossProduct(&Va, &Vb, &Vc); // 向量积: Va x Vb -> Vc
    end{verbatim}

    注意 Va, Vb, Vc 在向量积中向量元素个数须相同.


  • 单矩阵操作:
    CvMat *Ma, *Mb;
    
    cvTranspose(Ma, Mb); // transpose(Ma) -> Mb (不能对自身进行转置)
    CvScalar t = cvTrace(Ma); // trace(Ma) -> t.val[0]
    double d = cvDet(Ma); // det(Ma) -> d
    cvInvert(Ma, Mb); // inv(Ma) -> Mb

  • 非齐次线性系统求解:
    CvMat* A   = cvCreateMat(3,3,CV_32FC1);
    
    CvMat* x = cvCreateMat(3,1,CV_32FC1);
    CvMat* b = cvCreateMat(3,1,CV_32FC1);
    cvSolve(&A, &b, &x); // solve (Ax=b) for x

  • 特征值分析(针对对称矩阵):
    CvMat* A   = cvCreateMat(3,3,CV_32FC1);
    
    CvMat* E = cvCreateMat(3,3,CV_32FC1);
    CvMat* l = cvCreateMat(3,1,CV_32FC1);
    cvEigenVV(&A, &E, &l); // l = A的特征值 (降序排列)
    // E = 对应的特征向量 (每行)

  • 奇异值分解SVD:
    CvMat* A   = cvCreateMat(3,3,CV_32FC1);
    
    CvMat* U = cvCreateMat(3,3,CV_32FC1);
    CvMat* D = cvCreateMat(3,3,CV_32FC1);
    CvMat* V = cvCreateMat(3,3,CV_32FC1);
    cvSVD(A, D, U, V, CV_SVD_U_T|CV_SVD_V_T); // A = U D V^T

    标号使得 U 和 V 返回时被转置(若没有转置标号,则有问题不成功!!!).

视频序列操作

从视频序列中抓取一帧

  • OpenCV支持从摄像头或视频文件(AVI)中抓取图像.

  • 从摄像头获取初始化:
    CvCapture* capture = cvCaptureFromCAM(0); // capture from video device #0
    

  • 从视频文件获取初始化:
    CvCapture* capture = cvCaptureFromAVI("infile.avi");
    

  • 抓取帧:
    IplImage* img = 0;
    
    if(!cvGrabFrame(capture)){ // 抓取一帧
    printf("Could not grab a frame\n\7");
    exit(0);
    }
    img=cvRetrieveFrame(capture); // 恢复获取的帧图像

    要从多个摄像头同时获取图像, 首先从每个摄像头抓取一帧. 在抓取动作都结束后再恢复帧图像.

  • 释放抓取源:
    cvReleaseCapture(&capture);
    

    注意由设备抓取的图像是由capture函数自动分配和释放的. 不要试图自己释放它.

获取/设定帧信息

  • 获取设备特性:
    cvQueryFrame(capture); // this call is necessary to get correct
    
    // capture properties
    int frameH = (int) cvGetCaptureProperty(capture, CV_CAP_PROP_FRAME_HEIGHT);
    int frameW = (int) cvGetCaptureProperty(capture, CV_CAP_PROP_FRAME_WIDTH);
    int fps = (int) cvGetCaptureProperty(capture, CV_CAP_PROP_FPS);
    int numFrames = (int) cvGetCaptureProperty(capture, CV_CAP_PROP_FRAME_COUNT);

    所有帧数似乎只与视频文件有关. 用摄像头时不对,奇怪!!!.

  • 获取帧信息:
    float posMsec    =        cvGetCaptureProperty(capture, CV_CAP_PROP_POS_MSEC);
    
    int posFrames = (int) cvGetCaptureProperty(capture, CV_CAP_PROP_POS_FRAMES);
    float posRatio = cvGetCaptureProperty(capture, CV_CAP_PROP_POS_AVI_RATIO);

    获取所抓取帧在视频序列中的位置, 从首帧开始按[毫秒]算. 或者从首帧开始从0标号, 获取所抓取帧的标号. 或者取相对位置,首帧为0,末帧为1, 只对视频文件有效.

  • 设定所抓取的第一帧标号:
    // 从视频文件相对位置0.9处开始抓取
    
    cvSetCaptureProperty(capture, CV_CAP_PROP_POS_AVI_RATIO, (double)0.9);

    只对从视频文件抓取有效. 不过似乎也不成功!!!

存储视频文件

  • 初始化视频存储器:
    CvVideoWriter *writer = 0;
    
    int isColor = 1;
    int fps = 25; // or 30
    int frameW = 640; // 744 for firewire cameras
    int frameH = 480; // 480 for firewire cameras
    writer=cvCreateVideoWriter("out.avi",CV_FOURCC('P','I','M','1'),
    fps,cvSize(frameW,frameH),isColor);

    其他有效编码:

    CV_FOURCC('P','I','M','1')     = MPEG-1 codec
    
    CV_FOURCC('M','J','P','G') = motion-jpeg codec (does not work well)
    CV_FOURCC('M', 'P', '4', '2') = MPEG-4.2 codec
    CV_FOURCC('D', 'I', 'V', '3') = MPEG-4.3 codec
    CV_FOURCC('D', 'I', 'V', 'X') = MPEG-4 codec
    CV_FOURCC('U', '2', '6', '3') = H263 codec
    CV_FOURCC('I', '2', '6', '3') = H263I codec
    CV_FOURCC('F', 'L', 'V', '1') = FLV1 codec

    若把视频编码设为-1则将打开一个编码选择窗口(windows系统下).

  • 存储视频文件:
    IplImage* img = 0;
    
    int nFrames = 50;
    for(i=0;i cvGrabFrame(capture); // 抓取帧
    img=cvRetrieveFrame(capture); // 恢复图像
    cvWriteFrame(writer,img); // 将帧添加入视频文件
    }

    若想在抓取中查看抓取图像, 可在循环中加入下列代码:

    cvShowImage("mainWin", img);
    
    key=cvWaitKey(20); // wait 20 ms

    若没有20[毫秒]延迟,将无法正确显示视频序列.

  • 释放视频存储器:
    cvReleaseVideoWriter(&writer);

2008年7月19日星期六

冷笑话100则!

1:从前有个人钓鱼,钓到了只鱿鱼。
鱿鱼求他:你放了我吧,别把我烤来吃啊。
那个人说:好的,那么我来考问你几个问题吧。
鱿鱼很开心说:你考吧你考吧!
然后这人就把鱿鱼给烤了..
2:我曾经得过精神分裂症,但现在我们已经康复了。
3:一留学生在美国考驾照,前方路标提示左转,他不是很确定,问考官:
“turn left?”
答:“right”
于是……挂了..
4:有一天绿豆自杀从5楼跳下来,流了很多血,变成了红豆;一直流脓,又变成了黄豆;伤口结了疤,最后成了黑豆。
5:小明理了头发,第二天来到学校,同学们看到他的新发型,笑道:小明,你的头型好像个风筝哦!小明觉得很委屈,就跑到外面哭。哭着哭着~他就飞起来了…………
6:有个人长的像洋葱,走着走着就哭了…….
7:小企鹅有一天问他奶奶,“奶奶奶奶,我是不是一只企鹅啊?”“是啊,你当然是企鹅。”小企鹅又问爸爸,“爸爸爸爸,我是不是一只企鹅啊?”“是啊,你是企鹅啊,怎么了?”“可是,可是我怎么觉得那么冷呢?”
8:有一对玉米相爱了…
于是它们决定结婚…
结婚那天…
一个玉米找不到另一个玉米了…
这个玉米就问身旁的爆米花:你看到我们家玉米了吗?
爆米花:亲爱的,人家穿婚纱了嘛…….
9:音乐课上 老师弹了一首贝多芬的曲子
小明问小华:“你懂音乐吗?”
小华:“是的”
小明:“那你知道老师在弹什麼吗?”
小华: “钢琴。”
10:Q:有两个人掉到陷阱里了,死的人叫死人,活人叫什么?
A:叫救命啦!
11:提问:布和纸怕什么?
回答:布怕一万,纸怕万一。
原因:不(布)怕一万,只(纸)怕万一。
12:有一天有个婆婆坐车…
坐到中途婆婆不认识路了….
婆婆用棍子打司机屁股说:这是哪?
司机:这是我的屁股…..
13: 一个鸡蛋去茶馆喝茶,结果它变成了茶叶蛋;一个鸡蛋跑去松花江游泳,结果它变成了松花蛋;一有个鸡蛋跑到了山东,结果变成了鲁(卤)蛋;一个鸡蛋无家可 归,结果它变成了野鸡蛋;一个鸡蛋在路上不小心摔了一跤,倒在地上,结果变成了导弹;一个鸡蛋跑到人家院子里去了,结果变成了原子弹;一个鸡蛋跑到青藏高 原,结果变成了氢弹;一个鸡蛋生病了,结果变成了坏蛋;一个鸡蛋嫁人了,结果变成了混蛋;一个鸡蛋跑到河里游泳,结果变成了核弹;一个鸡蛋跑到花丛中去 了,结果变成了花旦;一个鸡蛋骑着一匹马,拿着一把刀,原来他是刀马旦;一个鸡蛋是母的,长的很丑,结果就变成了恐龙蛋;一个鸡蛋是公的,他老婆在外面和 别的鸡蛋通奸,结果他变成了王八蛋;一个鸡蛋……
14:主持人问:猫是否会爬树?老鹰抢答:会!主持人:举例说明!老鹰含泪:那年,我睡熟了,猫爬上了树…后来就有了猫头鹰…
15:俩屎壳螂讨论福利彩票,甲说:我要中了大奖就把方圆50里的厕所都买下来,每天吃个够!乙说:你丫太俗了!我要是中了大奖就包一活人,每天吃新鲜的!
16:why the chicken cross the street
答案 to get another side
17:甲:那个人在干什么?
乙:他在发抖。
甲:他为什么要发抖呢?
乙:他冷呀。
甲:哦,原来发抖就不会冷拉。
甲:……
18:有个香蕉先生和女朋友约会,走在街上,天气很热,香蕉先生就把衣服脱掉了,之后他的女朋友就摔倒了………
19:一个香肠被关在冰箱里
感觉很冷,然后看了看身边的另一根,有了点安慰,说:“看你都冻成这样了,全身都是冰!”结果那根说:“对不起,我是冰棒。”
20:.从前有一个棉花糖去打了球打了很长时间.他说:好累啊,我觉得我整个人都软下来了……….
21:这位跳水运动员的动作难度很大,他做了一个转体三周接前空翻三周半接后空翻一个月。
22:MM找大学迷路了。遇见一位文质彬彬的教授。
MM:请问,我怎样才能到大学去?
教授:只有努力读书,才可以上大学。
23:局长与科长共乘电梯,局长放一屁后对科长说:你放屁了!科长说:不是我放的…不久科长被免职,局长在会上说:屁大的事你都担待不起,要你何用?
24:小姐:现在生意不好做呀!
老大:为什么?
小姐:“禽流感…..”
25:一女遇劫匪颤抖曰:“俺是XX学校的,刚毕业,工作都没找到,真的没有钱……”
劫匪听后竟然痛哭流涕,“妹子,俺也是XX学校的,你拿好学生证,前面抢劫的还是XX学校的,你放心,阿拉绝不抢自己人!”
26:想和女友ML,女友曰不洗澡不行,应允天冷可洗“局部”,洗毕,女友极为娇羞道:“亲爱的,你好好懒呦,用哪洗哪……”偶听完晕倒,偶就是刷了个牙啊~~~(巨隐讳的冷笑话)
27:一个盲人乞丐戴着墨镜在街上行乞。
一个醉汉走过来,觉得他可怜,就扔了一百元给他。
走了一段路,醉汉一回头,恰好看见那个盲人正对着太阳分辨那张百元大抄的真假。
醉汉过来一把夺回钱道:“你TMD不想活了,竟敢骗老子!”
盲人乞丐一脸委屈说:“大哥,真对不起啊,我是替一个朋友在这看一下,他是个瞎子,去上厕所了,其实我是个哑巴。”
“哦,是这样子啊,”于是醉汉扔下钱,又摇摇晃晃地走了……
28:禽流感——都是“天屎”惹的祸!!!
有两种人得禽流感的几率极大——1.“禽兽” ;2.“禽兽不如”的人 …….
29:A:哎,你怎么学会抽烟了?
B:我从亚当夏娃偷吃禁果的时候就会了~
C:知道亚当夏娃为什么会偷吃禁果吗?
AB:不知!
C:因为亚当没有烟!(提示:谐音一个字)
30:某人刚被女友抛弃,碰巧在大街上撞见前女友和新欢调情,他越看越气,想羞辱他们一下。于是很有礼貌上前打了个招呼,并很鄙视地对女友新欢说:“我用过的旧货你也不嫌弃!”正当他为自己创意得意的时候,前女友却笑出声道:“外面一寸是旧的,里面全是崭新的!”
31:分手时,她给了我一个吻,那感觉——就好像人民日报一样真实……
32:刚刚看师姐的电脑屏幕上方有个类似新闻滚动条的东西,上面的文字过得非常快。
偶好奇问:这是歌词吗?
师姐:是呀!
师姐:怎么过得这么快?都没看清!
师姐:周杰伦的!!
33:妻:我真是瞎了眼踩到狗屎才会嫁给你。
夫:我才真是瞎了眼踩到狗屎才会娶你。
狗屎:我好倒霉喔!躺在那里都被你们俩给踩到……
34:高考化学题:A和B可以相互转化,B在沸水中可以生成C,C在空气 中氧化成D,D有臭鸡蛋气味,问A,B,C,D各是什么?
我答:A是鸡,B是生鸡蛋,C是熟鸡蛋,D当然是臭鸡蛋啦!
35:橡皮、老虎皮、狮子皮哪一个最不好?
答:橡皮。
因为橡皮擦(橡皮差)。
36:问:3个头一只脚的是什么东西???
答案:3个头一只脚的怪物!!!!!!
37:蚂蚁去沙漠,为什么沙子上没有留下他的脚印,而只留下一条线呢?
答案:因为它是骑脚踏车的!
蚂蚁从沙漠回家了,他没有通知任何人,但是他家人却知道他回来了!为什么啊!
答案:看见他停在楼下的脚踏车…….
38:有一天一个女吸毒犯被抓到警局,police看见她的手上有刺青,就问她你干嘛把你男朋友的名字刺在手上,他叫小良是不是…啊..是不是.快说,说..他有没有吸毒阿….快说
只见那个女吸毒犯抬起头带着愤怒的眼神
对police说
这是恨啦….
40:一天,小美和她男友开车出去兜风,
车快没油了,刚好旁边有个加油站,开过去的时候,突然一阵狂风把她男友的帽子刮跑了。
小美的男友对她说:
「我去捡帽子,你帮我加油。」
男友刚跑开不远,就听到小美在他后面大喊:
「加油!加油!」
41:一只猩猩经过树林,不小心采到了长臂猿的粪便,
好心的猩猩把猿打扫了分辩.
过了不久他们相爱了,别人问你们是怎么走到一起的?
猩猩回答说:”是猿粪(缘分).!”
42::有一个胖子……….
从高楼跳下…
结果变成了…….
死胖子..
43:有一只鸭子叫小黄,有一天它过马路时被车撞了一下,大叫:“呱!”从此它就变成了小黄瓜……
44: 有一只企鹅,他的家离北极熊家特别远,要是靠走的话,得走20年才能到。有一天,企鹅在家里呆着特别无聊,准备去找北极熊玩,与是他出门了,可是走到路的 一半的时候发现自己忘记锁门了,这就已经走了10年了,可是门还是得锁啊,于是企鹅又走回家去锁门。锁了门以后,企鹅再次出发去找北极熊,等于他花了40 年才到了北极熊他们家……然后企鹅就敲门说:“北极熊北极熊,企鹅找你玩来了!”结果北极熊开门以后你猜他说什么?“还是去你家玩吧~”
45: 小白兔蹦蹦跳跳到面包房,问:“老板,你们有没有一百个小面包啊?” 老板:“啊,真抱歉,没有那么多”“这样啊。。。”小白兔垂头丧气地走了。第二天,小白兔蹦蹦跳跳到面包房,“老板,有没有一百个小面包啊?” 老板:“对不起,还是没有啊”“这样啊。。。”小白兔又垂头丧气地走了。第三天,小白兔蹦蹦跳跳到面包房,“老板,有没有一百个小面包啊?”老板高兴的 说:“有了,有了,今天我们有一百个小面包了!!” 小白兔掏出钱:“太好了,我买两个!”
46:小明说:「阿康,问你“有一只鲨鱼吃下了一颗绿豆,结果它变成了什么“?」 阿康说:「我不知道,答案是什么?」小明说:「嘿!嘿!答案是“绿豆沙(绿豆鲨)“,你很笨喔!」
47:老师问一同学怎么减少白色污染? 同学答:把饭盒做成蓝色.
48:有个人,他肠胃不好.一天,他来到胃病医院看病,对医生说:“我吃什么拉什么,吃西 瓜拉西瓜,吃黄瓜拉黄瓜!“医生想了想,对他说:“我看你只有吃屎了!”
49:飞机上,一位空中小姐问一个小女孩说:“为什么飞机飞这么高都不会撞到星星呢?“ 小女孩回答到:“我知道,因为星星会’闪’啊!”
50:有一只北极熊和一只企鹅在一起耍,企鹅把身上的毛一根一根地拔了下来,拔完之后,对北极熊说:“好冷哦!“北极熊听了,也把自己身上的毛一根一根地拔了下来,转头对企鹅说:“果然很冷!”
51:Q:非洲食人族的酋长吃什么?
A:人啊!
Q:那有一天,酋长病了,医生告诉他要吃素,那他吃什么?
A:吃植物人!
52:冰箱里有两根香肠,过了很久,
一跟香肠抖了一下,哇!好冷啊~!
另一根香肠十分惊奇地说,咦?你是香肠怎么会说话?
53:有一天,
有只公鹿越跑越快,
跑到最后 ,
它就变成高速公鹿了。
54:有一天,老师带一群小朋友到山上采水果,
她宣布说:“小朋友,采完水果后,我们统一一起洗,洗完可以一起吃。”
所有小朋友都跑去采水果了。
集合时间一到,所有小朋友都集合了。
老师:“小华,你采到什么?”
小华:“我在洗苹果,因为我采到苹果。”
老师:“小美你呢?”
小美:“我在洗蕃茄,因为我采到蕃茄。”
老师:“小朋友都很棒哦!那阿明你呢?”
阿明:“我在洗布鞋,因为我踩到大便。”
55:老师在课堂上对小明提问,小明站起来却一声不吭。
老师:小明?
老师:小明??
老师:小明!你怎么回事啊?你到底知不知道答案啊?好歹吱一声啊!
小明:吱~
56:一只大象问骆驼:‘你的咪咪怎么长在背上?’
骆驼说:‘死远点,我不和鸡鸡长在脸上的东西讲话!
57:如何让饮料变大杯?
念大悲咒
58:小明:今天几度阿?
小华:零下3度阿!
小明:难怪这麼冷。
59:有个小男孩放学回家从窗外窥见一个女人躺在床上狂揉胸部喊到我要男人我要男人!
第二天小男孩再从窗外过时发现女人身上躺了个男人,
于是小男孩回家躺在床上狂揉胸部喊到我要自行车我要自行车!
60:从前从前有一只鸟,
他每天都会经过一片玉米田,
但是很不幸的,
有一天那片玉米田发生了火灾,
所有的玉米都变成了爆米花!!!
小鸟飞过去以后……
以为下雪,就冷死了…
61:说有一只北极熊,因为雪地太刺眼了,必须要戴墨镜才能看东西,
可是他找不到墨镜,于是闭着眼睛爬来爬去在地上找,爬呀爬呀,把 手脚都爬的脏兮兮的才找到墨镜。戴上墨镜,对着镜子一照,这才发现:哦,原来我是一只熊猫。
62:自然课老师问:为什么人死后身体是冷的?
没人回答。
老师又问:没人知道吗?
这时,有个同学站起来说:那是因为心静自然凉。
63:小明在一次车祸中失去了一条腿,
小明在一次车祸中又失去了一条腿,
又一次车祸中小明失去了他的另一条腿,
一次车祸中小明又失去了他的一条腿,
其实小明是一条狗.
64:一天,A,B,C三个人一起出去玩,走在路上闲晃了很久。
后来A就说,好无聊,我好想去打B。
然后C看了A一眼, 就把B拖到巷子里去打。
65:三只小兔拉便便
第一只是长条的 。
第二只是圆球的。
第三只居然是三角形的 。
问,它答:我用手捏的。
66:台湾什么时候会想要统一?
买方便面的时候
67:有一天小强问他爸爸:“爸爸,我是不是傻孩子啊?”
爸爸说:“傻孩子,你怎么会是傻孩子呢?”
68:小明回家时,隔壁的狗突然跑出来咬他,他一气之下拿起竹子要打它,
狗的主人看到小明打他的狗,就不高兴的说:打狗也要看主人,没听过吗?
这时小明就说:好!我会一边看着你,一边打你家的狗。
69:虫虫:小花,你用我的铅笔了吗?
小花:没有,我没用。
虫虫:你真没用?
小花:我真没用!
虫虫:唉,你是第17个承认自己没用的人了
70:蚂蚁从喜玛拉雅山上摔下来后是怎么死的?
答案:饿死的。因为太轻~所以飘下来要很久…
80:为什么小狗越变越小?
答:因为它越走越远。
81:从前,有一只马!它跑着跑着就掉进海里。
所以,它变成了一只“海马”!
这只马的另外一只马朋友,为了要去找掉到海里的马,结果却掉到河里。后来,他就就变成“河马”。
第三只马是只白马。它为了要找失踪的两个朋友,来到了交通混乱的城市。
它连续被好几台车子给辗过,使得身上出现好几条黑条纹。
结果,它变成“斑马”了!
第四只马为了找寻前面三个的同伴,有一天,它来到一间工厂,结果被改造成“铁马”。
但后来,那些马还是难逃被吃的命运,通通被作成了“沙其马”,肆虐所及,所有马儿无一幸免,成了一个无马的世界……
然后,有一群人看到这篇笑话后忍不住的说:“马的~真冷”。
最后,为了纪念这个笑话,有人将它编订成课,我们叫它“马赛课”!
82:小明欠地下钱庄20万,小明苦苦哀求他多宽限几天,
钱庄的人说:明天一定要还,不然的话……,剁掉2只手指;
后天的话……,在剁4只;第3天的话……
小明:是不是不用还了
钱庄的人:NO,到时候你就变成小叮当了。
83:有个人一天碰到上帝
上帝突然大发善心打算给那人一个愿望
上帝问:你有什么愿望吗?
那个人想了想说: 听说猫都有9条命,那请您赐给我9条命吧!
上帝说:你的愿望实现咯!
一天,那个人闲着无聊,
想说去死一死算了, 反正有9条命嘛
就躺在铁轨上,
结果一辆火车开过去,
那人还是死了.
这是为什么呢?
因为那台火车的车厢有10节.
84:一个家伙到医院去检查,并做了许多测试。
医生说:有好消息、也有坏消息!看过你的测试结果后,我发现你有潜在的同性恋倾向!!而且难以根治!
这个家伙说:我的天啊!那好消息呢?
医生腼腆的说:我发现你还蛮可爱的耶
85:一个猎人带着猎狗去打猎,在林子里溜了一天都没有猎物。
天黑了,不甘心的他还是不停骑马在林子里转,
马忽然说:‘你都不让我休息,想累死我啊!?’
猎人听到吓了一跳,立刻从马背上滚下来,拉着猎狗就逃跑,跑到一课大树下喘气时,狗拍拍胸口对他说:‘吓死我了,马居然会说话!’
于是猎人当场被吓死了
86:狼、老虎和狮子谁玩游戏一定会被淘汰? 狼
因为:桃太郎(淘汰狼)
87:一天A拣了一面镜子对着镜子照了照说;这里边的人好面熟啊
B说;是吗?我看看(接过镜子),我啊!我你都不认识了啊?
88:番茄A和番茄B去逛街。
B问A:我们去哪?
A不回答。
B又问:我们去哪?
A还是不回答。
B又问了一次。
番茄A转过头对番茄B说:我们不是番茄来的吗?为什么我们会说话呢?
89:从前,有一只白猫和一只黑猫
一天
白猫掉到水里去了
黑猫把 它救了上来
白猫对黑猫说了一句话
“Q:这句话是什么?
.
.
.
.
.
.
.
.
.
.
.
.
..
“喵”
90:A:“你知道我昨天晚上在网吧干嘛吗?”
B:“在干吗;”
A:“上网呗;”
B:“。。。”
91:两只苍蝇去吃饭。
小的问大的:大哥,为什么我们每天都要吃屎?
大的说:吃饭的时候不要说这么恶心的东西!!
92:草船中
鲁肃:“这样真的可以借到箭吗?孔明先生?”
诸葛亮:“相信我。”
鲁肃:“可是我还是有些担心……”
诸葛亮:“没必要。”
鲁肃:“可是,你不觉得船里越来越热么?”
诸葛亮:“这么说起来是有一点碍…有什么不对劲吗?”
鲁肃:“是啊,我担心敌人射的是火箭……”
诸葛亮:“哎!?子敬 ̄ ̄你会游泳么 ̄ ̄ ̄我不会 ̄ ̄ ̄”
93:一猴子吃花生前都要先塞进屁股再拿出来吃。
对此管理员解释道:曾有人喂它桃 子,
结果桃核拉不出来,猴子吓怕了,现在一定要量好再吃。
94:医院为防止病人出逃外设100道,两精神病患者仍欲逃出医院。于夜黑努力
翻墙。至第30道墙下,
“累了么?” ,
“不累。”于是二人继续向外翻。
至第60道墙下,
“你累了么?”
“不累。”于是二人继续向外翻,
至第99道墙下,
“你累了么? ”
“累了”
“那好,我们翻回去吧”
95:小明:在某个溪边,有大宝、大雄、大志、大伟共四个男孩脱光光在玩水,
突然有人在溪边电鱼,这四个男孩都被电到了!猜一种电器用品。
阿康:嗯…. 不知道耶~
小明:答案是”电视机”(电四鸡)!嘿嘿!
96:小骆:爸爸,为什麼我们要有驼峰呢?
驼爸:因为沙漠中没有水,有驼峰才可以储存水分啊!
小骆:爸爸,为什麼我们要有长长的毛呢?
驼爸:因为沙漠中风沙大,我们必须靠它阻挡风砂,才看得见啊!
小骆:爸爸,为什麼我们要有厚厚的蹄呢?
驼爸:因为沙漠中都是沙,这样我们才站得稳啊!
小骆:爸爸,最后一个问题,那我们在动物园干嘛呢?
97:母鸡在孵蛋,有个蛋从它屁屁钻出来了
母鸡:“你干吗?”
鸡蛋:“你放屁好臭……”
98:有个人的名字叫“杜子藤”
老师点名时问
“杜子藤呢?”
同学说:“他肚子疼。”
99:我女朋友约我去她家看电影。到了她家之后,
她用签字笔在墙上写了‘电影’两个字,
我们两个就坐在马桶上看了起来。
100:某清晨,以严厉著称的某长官问晨练小兵:“你冷吗?”
小兵答:“不冷!”
长官恼:“那你颤颤什么?”
小兵答:“冻的!”

2008年7月5日星期六

KLT

KLT: An Implementation of the
Kanade-Lucas-Tomasi Feature Tracker


KLT is an implementation, in the C programming language, of a feature tracker for the computer vision community. The source code is in the public domain, available for both commercial and non-commerical use.

The tracker is based on the early work of Lucas and Kanade [1], was developed fully by Tomasi and Kanade [2], and was explained clearly in the paper by Shi and Tomasi [3]. Later, Tomasi proposed a slight modification which makes the computation symmetric with respect to the two images -- the resulting equation is derived in the unpublished note by myself [4]. Briefly, good features are located by examining the minimum eigenvalue of each 2 by 2 gradient matrix, and features are tracked using a Newton-Raphson method of minimizing the difference between the two windows. Multiresolution tracking allows for relatively large displacements between images. The affine computation that evaluates the consistency of features between non-consecutive frames [3] was implemented by Thorsten Thormaehlen several years after the original code and documentation were written.

Some Matlab interface routines: klt_read_featuretable.m

Note: An alternate Lucas-Kanade implementation can be found in Intel's OpenCV library. This implementation, described in the note by Bouguet, does a better job of handling features near the image borders, and it is more computationally efficient (approximately 30% on my desktop system). However, it does not contain the affine consistency check. Another alternative is GPU_KLT, which is an implementation of KLT for a graphics processing unit (GPU), which speeds up the run time considerably. A Matlab implementation of a single template tracker is available at Lucas-Kanade 20 Years On.

References

[1] Bruce D. Lucas and Takeo Kanade. An Iterative Image Registration Technique with an Application to Stereo Vision. International Joint Conference on Artificial Intelligence, pages 674-679, 1981.

[2] Carlo Tomasi and Takeo Kanade. Detection and Tracking of Point Features. Carnegie Mellon University Technical Report CMU-CS-91-132, April 1991.

[3] Jianbo Shi and Carlo Tomasi. Good Features to Track. IEEE Conference on Computer Vision and Pattern Recognition, pages 593-600, 1994.

[4] Stan Birchfield. Derivation of Kanade-Lucas-Tomasi Tracking Equation. Unpublished, January 1997.


Website maintained by Stan Birchfield

2008年7月3日星期四

如何用摄像头来测距(opencv)



原 理

假设激光束是与摄像头的光轴完全平行,激光束的中心落点在在摄像头的视域中是最亮的点。激光束照射到摄像头视域中的跟踪目标上,那么摄像头可以捕捉到这个 点,通过简单的图像处理的方法,可以在这侦图像中找到激光束照射形成的最亮点,同时可以计算出Y轴上方向上从落点到图像中心的象素的个数。这个落点越接近 图像的中心,被测物体距离机器人就越远。由下图图可以计算距离D:

(1)

等式中h是一个常量,是摄像头与激光发射器之间的垂直距离,可以直接测量获得。

θ可通过下式计算:
θ=Num*Rop+Offset (2)
其中:Num是从图像中心到落点的像素个数
Rop是每个像素的弧度值
Offset是弧度误差
合并以上等式可以得到:
(3)
Num可以从图像上计算得到。Rop和Offset需要通过实验计算获得。首先测量出D的准确值,然后根据等式(1)可以计算出准确的θ,根据 等式(2)可到只含有参数Rop和Offset的方程。在不同的距离多次测量D的准确值计算θ,求解方程组可以求出Rop和Offset。这里 Rop=0.0030354,Offset=0.056514344。


程 序
头文件:
class LaserRange
{
public:
struct RangeResult * GetRange(IplImage * imgRange,IplImage * imgDst);
LaserRange();
virtual ~LaserRange();
private:
unsigned int maxW;
unsigned int maxH;
unsigned int MaxPixel;
RangeResult * strctResult;

// Values used for calculating range from captured image data
const double gain;
// Gain Constant used for converting pixel offset to angle in radians
const double offset; // Offset Constant
const double h_cm; // Distance between center of camera and laser
unsigned int pixels_from_center;
// Brightest pixel location from center
void Preprocess(void * img,IplImage * imgTemp);
};
cpp文件:

LaserRange::LaserRange():gain(0.0030354),offset(0),h_cm(4.542)
{
maxW=0;
maxH=0;
MaxPixel=0;

pixels_from_center=0;
// Brightest pixel location from center
strctResult=new RangeResult;

strctResult->maxCol=0;
strctResult->maxRow=0;
strctResult->maxPixel=0;
strctResult->Range=0.0;
}
LaserRange::~LaserRange()
{
if(NULL!=strctResult) delete strctResult;
}
struct RangeResult * LaserRange::GetRange(IplImage * imgRange,IplImage * imgDst)
{
if(NULL==imgRange) return strctResult;
Preprocess(imgRange,imgDst);

pixels_from_center = abs(120-maxH);
// Calculate range in cm based on bright pixel location, and setup specific constants
strctResult->Range= h_cm/tan(pixels_from_center * gain + offset);

strctResult->PixfromCent=pixels_from_center;
strctResult->maxCol=maxW;
strctResult->maxRow=maxH;
strctResult->maxPixel=MaxPixel;
//strctResult->Range=0.0;
return strctResult;
}

void LaserRange::Preprocess(void *img, IplImage * imgTemp)
{
MaxPixel=0;
//处理下一帧前 最大像素值清零;
IplImage* image = reinterpret_cast(img);

cvCvtPixToPlane( image,0 ,0 ,imgTemp , 0);

for( int j=((imgTemp->width-60)/2-1); j<(imgTemp->width-40)/2+59; j++)
{
for(int i=5; i
height-5; i++)
{

if((imgTemp->imageData[i*imgTemp->widthStep+j])>MaxPixel)
{
if( ((imgTemp->imageData[(i-1)*imgTemp->widthStep+j])>MaxPixel) && ((imgTemp->imageData[(i-1)*imgTemp->widthStep+j+1])>MaxPixel) &&((imgTemp->imageData[(i-1)*imgTemp->widthStep+j-1])>MaxPixel) )
{
if( ((imgTemp->imageData[(i+1)*imgTemp->widthStep+j])>MaxPixel) && ((imgTemp->imageData[(i+1)*imgTemp->widthStep+j+1])>MaxPixel) &&((imgTemp->imageData[(i+1)*imgTemp->widthStep+j-1])>MaxPixel) )
{
if((imgTemp->imageData[i*(imgTemp->widthStep)+j+1])>MaxPixel)
{
if((imgTemp->imageData[i*(imgTemp->widthStep)+j-1])>MaxPixel)
{
MaxPixel=imgTemp->imageData[i*imgTemp->widthStep+j] ;
maxW=j;
maxH=i;
}
}
}
}
}
}

}
调用函数:
int CLaserVisionDlg::CaptureImage()
{
// CvCapture* capture = 0;

// capture = cvCaptureFromCAM(0); //0表示设备号
if( !capture )
{
fprintf(stderr,"Could not initialize capturing...\n");
return -1;
}

// cvNamedWindow( "LaserRangeImage", 1 );
// cvvNamedWindow( "image", 1);
cvvNamedWindow( "Dimage", 1);

for(;;)
{
IplImage* frame = 0;

if(isStop) break;
frame = cvQueryFrame( capture ); //从摄像头抓取一副图像框架
if( !frame )
break;
if( !imgOrign )
{
//allocate all the buffers
imgOrign = cvCreateImage( cvGetSize(frame), 8, 3 ); //创建一副图像
imgOrign->origin = frame->origin;

}
cvCopy( frame, imgOrign, 0 ); //将图frame复制到image
//cvShowImage("LaserRangeImage",imgOrign);


if(!imgDest)
{
imgDest=cvCreateImage( cvSize( imgOrign->width,imgOrign->height),8,1);
cvZero( imgDest );
}
struct RangeResult * temp= laservsion.GetRange(imgOrign,imgDest);
cvLine( imgOrign,cvPoint(temp->maxCol,0), cvPoint(temp->maxCol,imgOrign->height),cvScalar(100,100,255,0),1,8,0);
cvLine( imgOrign,cvPoint(0,temp->maxRow), cvPoint(imgOrign->width,temp->maxRow),cvScalar(100,100,255,0),1,8,0);


// cvvShowImage( "image", imgOrign);
cvSaveImage("image.bmp", imgOrign);

cvvShowImage( "Dimage", imgDest);

//在PictureBox上显示图片
CDC* pDC = GetDlgItem(IDC_Picture)->GetDC();
CDC dcmemory;
BITMAP bm;
dcmemory.CreateCompatibleDC(pDC);
CBitmap* pBmp;
CString szFileName = "image.bmp";
HBITMAP hBk = (HBITMAP)::LoadImage(NULL,szFileName,IMAGE_BITMAP,0,0,LR_LOADFROMFILE);
if(NULL!=hBk)
{
pBmp=CBitmap::FromHandle(hBk);
pBmp->GetObject(sizeof(BITMAP), &bm);
dcmemory.SelectObject(pBmp);
pDC->BitBlt(0, 0, bm.bmWidth, bm.bmHeight, &dcmemory, 0, 0, SRCCOPY);
}


char str[80];
// To print message
CDC *pDCp= GetDC();
char str2[80];

// Display frame coordinates as well as calculated range
sprintf(str, "Pix Max Value=%d At x= %u, y= %u, PixfromCent= %d",temp->maxPixel,temp->maxCol, temp->maxRow, temp->PixfromCent);
sprintf(str2, "Range= %f cm ",temp->Range);
pDCp->TextOut(30, 33, str);
pDCp->TextOut(50, 50, str2);
ReleaseDC(pDCp);

int c = cvWaitKey(10);
// if( c == 'q' )
// break;
}
//cvReleaseCapture( &capture );
//cvDestroyWindow("LaserRangeImage");
// cvDestroyWindow( "image");
cvDestroyWindow( "Dimage");
return 0;
}

2008年7月1日星期二

利用光流法计算人体运动的速度与方向

1.方向的计算
首先计算图像各个象素的光流(opencv LK),然后建立4*4窗口对X,Y方向分别做统计求和,
然后求得 atan(yy/xx)作为光流方向,即为运动的方向.
2.速度的计算
利用帧差分得到运动图像,然后建立4*4窗口对图像进行统计求和,求和值作为权重,表示速度的比例.
即运动区域白色(255)面积越大,速度越大.

3.结果
大部分运动方向计算正确,少部分有错误,还需要改进算法.(利用统计?)

4.代码:

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 ;yheight; y++)
for(int x =0;xwidth;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)
alpha_angle = pi/2;
else
alpha_angle = abs(atan(yy/xx));

if(xx<0>0) alpha_angle = pi - alpha_angle ;
if(xx<0 alpha_angle =" pi">0 && yy<0) alpha_angle =" 2*pi">



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(x + scale_factor*cos(alpha_angle));
pt2.y = static_cast(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;

}

模糊数学(Fuzzy mathematics)及其应用

模糊数学(Fuzzy mathematics)及其应用
一 :引言
有一个古老的希腊悖论,是这样说的:“一粒种子肯定不叫一堆,两粒也不是,三粒也不是……另一方面,所有的人都同意,一亿粒种子肯定叫一堆。那么,适当的界限在哪里?我们能不能说,123585粒种子不叫一堆而123586粒就构成一堆?”
确 实,“一粒”和“一堆”是有区别的两个概念。但是,它们的区别是逐渐的,而不是突变的,两者之间并不存在明确的界限。换句话说,“一堆”这个概念带有某种 程度的模糊性。类似的概念,如“年老”、“高个子”、“年轻人”、“很大”、“聪明”、“漂亮的人”、“价廉物美”等等,不胜枚举。
经 典集合论中,在确定一个元素是否属于某集合时,只能有两种回答:“是”或者“不是”。我们可以用两个值0或1加以描述,属于集合的元素用1表示,不属于集 合的元素用0表示。然而上面提到的“年老”、“高个子”、“年轻人”、“很大”、“聪明”、“漂亮的人”、“价廉物美” 等情况要复杂得多。假如规定身高1.8米算属于高个子范围,那么,1.79米的算不算?照经典集合论的观点看:不算。但这似乎很有些悖于情理。如果用一个 圆,以圆内和圆周上的点表示集A,而且圆外的点表示不属于A。A的边界显然是圆周。这是经典集合的图示。现在,设想将高个子的集合用图表示,则它的边界将 是模糊的,即可变的。因为一个元素(例如身高1.75米的人)虽然不是100%的高个子,却还算比较高,在某种程度上属于高个子集合。这时一个元素是否属 于集合,不能光用0和1两个数字表示,而可以取0和1之间的任何实数。例如对1.75米的身高,可以说具有70%属于高个子集合的程度。这样做似乎罗嗦, 但却比较合乎实际。
精 确和模糊,是一对矛盾。根据不同情况有时要求精确,有时要求模糊。比如打仗,指挥员下达命令:“拂晓发起总攻。”这就乱套了。这时,一定要求精确:“×月 ×日清晨六时正发起总攻。”我们在一些旧电影中还能看到各个阵地的指挥员在接受命令前对对表的镜头,生怕出个半分十秒的误差。但是,物极必反。如果事事要 求精确,人们就简直无法顺利的交流思想——两人见面,问:“你好吗?”可是,什么叫“好”,又有谁能给“好”下个精确的定义?
有些现象本质上就是模糊的,如果硬要使之精确,自然难以符合实际。例如,考核学生成绩,规定满60分为合格。但是,59分和60分之间究竟有多大差异,仅据1分之差来区别及格和不及格,其根据是不充分的。
不 仅普遍存在着边界模糊的集合,就是人类的思维,也带有模糊的特色。有些现象是精确的,但是,适当的模糊化可能使问题得到简化,灵活性大为提高。例如,在地 里摘玉米,若要找一个最大的,那很麻烦,而且近乎迂腐。我们必须把玉米地里所有的玉米都测量一下,再加以比较才能确定。它的工作量跟玉米地面积成正比。土 地面积越大,工作越困难。然而,只要稍为改变一下问题的提法:不要求找最大的玉米,而是找比较大的,即按通常的说法,到地里摘个大玉米。这时,问题从精确 变成了模糊,但同时也从不必要的复杂变成意外的简单,挑不多的几个就可以满足要求。工作量甚至跟土地无关。因此,过分的精确实际成了迂腐,适当的模糊反而 灵活。
显 然,玉米的大小,取决于它的长度、体积和重量 。大小虽是模糊概念,但长度、体积、重量等在理论上都可以是精确的。然而,人们在实际判断玉米大小时,通常并不需要测定这些精确值。同样,模糊的“堆”的 概念是建立在精确的“粒”的基础上,而人们在判断眼前的东西叫不叫一堆时,从来不用去数“粒”。有时,人们把模糊性看成一种物理现象。近的东西看得清,远 的东西看不清,一般的说,越远越模糊。但是,也有例外的情况:站在海边,海岸线是模糊的;从高空向下眺望,海岸线却显得十分清晰。太高了,又模糊。精确与 模糊,有本质区别,但又有内在联系,两者相互矛盾、相互依存也可相互转化。所以,精确性的另一半是模糊。
对 模糊性的讨论,可以追溯得很早。20世纪的大哲学家罗素(B.Russel)在1923年一篇题为《含糊性》(Vagueness)的论文里专门论述过我 们今天称之为“模糊性”的问题(严格地说,两者稍有区别),并且明确指出:“认为模糊知识必定是靠不住的,这种看法是大错特错的。”尽管罗素声名显赫,但 这篇发表在南半球哲学杂志的文章并未引起当时学术界对模糊性或含糊性的很大兴趣。这并非是问题不重要,也不是因为文章写得不深刻,而是“时候未到”。罗素 精辟的观点是超前的。长期以来,人们一直把模糊看成贬义词,只对精密与严格充满敬意。20世纪初期社会的发展,特别是科学技术的发展,还未对模糊性的研究 有所要求。事实上,模糊性理论是电子计算机时代的产物。正是这种十分精密的机器的发明与广泛应用,使人们更深刻地理解了精密性的局限,促进了人们对其对立 面或者说它的“另一半”——模糊性的研究。
集 合是现代数学的基础,模糊集合一提出,“模糊”观念也渗透到许多数学分支。模糊数学的发展速度也是相当快的。从发表的论文看,几乎是指数般的增长。模糊数 学的研究可分三个方面:一是研究模糊数学的理论,以及它和精确数学、统计数学的关系;二是研究模糊语言和模糊逻辑;三是研究模糊数学的应用。在模糊数学的 研究中,目前已有模糊拓扑学、模糊群论、模糊凸论、模糊概率、模糊环论等分支。虽然模糊数学是一门新兴学科,但它已初步应用于自动控制、模式识别、系统理 论、信系检索、社会科学、心理学、医学和生物学等方面。将来还可能出现模糊逻辑电路、模糊硬件、模糊软件和模糊固件,出现能和人用自然语言对话、更接近于 人的智能的新的一类计算机。所以,模糊数学将越来越显示出它的巨大生命力。
二: 模糊数学的产生
二 十世纪六十年代,产生了模糊数学这门新兴学科。 现代数学是建立在集合论的基础上。集合论的重要意义就一个侧面看,在与它把数学的抽象能力延伸到人类认识 过程的深处。一组对象确定一组属性,人们可以通过说明属性来说明概念(内涵),也可以通过指明对象来说明它。符合概念的那些对象的全体叫做这个概念的外 延,外延其实就是集合。从这个意义上讲,集合可以表现概念,而集合论中的关系和运算又可以表现判断和推理,一切现实的理论系统都一可能纳入集合描述的数学 框架。
但是,数学的发展也是阶段性的。经典集合论只能把自己的表现力限制在那些有明确外延的概念和事物上,它明确地限定:每个集合都必须 由明确的元素构成,元素对集合的隶属关系必须是明确的,决不能模棱两可。对于那些外延不分明的概念和事物,经典集合论是暂时不去反映的,属于待发展的范 畴。
在较长时间里,精确数学及随机数学在描述自然界多种事物的运动规律中,获得显著效果。但是,在客观世界中还普遍存在着大量的模糊现象。以前人们回避它,但是,由于现代科技所面对的系统日益复杂,模糊性总是伴随着复杂性出现。
各门学科,尤其是人文、社会学科及其它“软科学”的数学化、定量化趋向把模糊性的数学处理问题推向中心地位。更重要的是,随着电子计算机、控制论、系统科学的迅速发展,要使计算机能像人脑那样对复杂事物具有识别能力,就必须研究和处理模糊性。
我们研究人类系统的行为,或者处理可与人类系统行为相比拟的复杂系统,如航天系统、人脑系统、社会系统等,参数和变量甚多,各种因素相互交错,系统很复杂,它的模糊性也很明显。从认识方面说,模糊性是指概念外延的不确定性,从而造成判断的不确定性。
在 日常生活中,经常遇到许多模糊事物,没有分明的数量界限,要使用一些模糊的词句来形容、描述。比如,比较年轻、高个、大胖子、好、漂亮、善、热、远……。 在人们的工作经验中,往往也有许多模糊的东西。例如,要确定一炉钢水是否已经炼好,除了要知道钢水的温度、成分比例和冶炼时间等精确信息外,还需要参考钢 水颜色、沸腾情况等模糊信息。因此,除了很早就有涉及误差的计算数学之外,还需要模糊数学。
人与计算机相比,一般来说,人脑具有处理模糊 信息的能力,善于判断和处理模糊现象。但计算机对模糊现象识别能力较差,为了提高计算机识别模糊现象的能力,就需要把人们常用的模糊语言设计成机器能接受 的指令和程序,以便机器能像人脑那样简洁灵活的做出相应的判断,从而提高自动识别和控制模糊现象的效率。这样,就需要寻找一种描述和加工模糊信息的数学工 具,这就推动数学家深入研究模糊数学。所以,模糊数学的产生是有其科学技术与数学发展的必然性。
模糊数学的研究内容
1965年,美国控制论专家、数学家查德发表了论文《模糊集合》,标志着模糊数学这门学科的诞生。
模糊数学的研究内容主要有以下三个方面:
第 一,研究模糊数学的理论,以及它和精确数学、随机数学的关系。察德以精确数学集合论为基础,并考虑到对数学的集合概念进行修改和推广。他提出用“模糊集合 ”作为表现模糊事物的数学模型。并在“模糊集合”上逐步建立运算、变换规律,开展有关的理论研究,就有可能构造出研究现实世界中的大量模糊的数学基础,能 够对看来相当复杂的模糊系统进行定量的描述和处理的数学方法。
在模糊集合中,给定范围内元素对它的隶属关系不一定只有“是”或“否”两种 情况,而是用介于0和1之间的实数来表示隶属程度,还存在中间过渡状态。比如“老人”是个模糊概念,70岁的肯定属于老人,它的从属程度是 1,40岁的人肯定不算老人,它的从属程度为 0,按照查德给出的公式,55岁属于“老”的程度为0.5,即“半老”,60岁属于“老”的程度0.8。查德认为,指明各个元素的隶属集合,就等于指定了 一个集合。当隶属于0和1之间值时,就是模糊集合。
第二,研究模糊语言学和模糊逻辑。人类自然语言具有模糊性,人们经常接受模糊语言与模糊信息,并能做出正确的识别和判断。
为了实现用自然语言跟计算机进行直接对话,就必须把人类的语言和思维过程提炼成数学模型,才能给计算机输入指令,建立和是的模糊数学模型,这是运用数学方法的关键。查德采用模糊集合理论来建立模糊语言的数学模型,使人类语言数量化、形式化。
如 果我们把合乎语法的标准句子的从属函数值定为1,那么,其他文法稍有错误,但尚能表达相仿的思想的句子,就可以用以0到1之间的连续数来表征它从属于“正 确句子”的隶属程度。这样,就把模糊语言进行定量描述,并定出一套运算、变换规则。目前,模糊语言还很不成熟,语言学家正在深入研究。
人们的思维活动常常要求概念的确定性和精确性,采用形式逻辑的排中律,既非真既假,然后进行判断和推理,得出结论。现有的计算机都是建立在二值逻辑基础上的,它在处理客观事物的确定性方面,发挥了巨大的作用,但是却不具备处理事物和概念的不确定性或模糊性的能力。
为了使计算机能够模拟人脑高级智能的特点,就必须把计算机转到多值逻辑基础上,研究模糊逻辑。目前,模糊罗基还很不成熟,尚需继续研究。
第 三,研究模糊数学的应用。模糊数学是以不确定性的事物为其研究对象的。模糊集合的出现是数学适应描述复杂事物的需要,查德的功绩在于用模糊集合的理论找到 解决模糊性对象加以确切化,从而使研究确定性对象的数学与不确定性对象的数学沟通起来,过去精确数学、随机数学描述感到不足之处,就能得到弥补。在模糊数 学中,目前已有模糊拓扑学、模糊群论、模糊图论、模糊概率、模糊语言学、模糊逻辑学等分支。
模糊数学的应用
模糊数学是一门新兴学 科,它已初步应用于模糊控制、模糊识别、模糊聚类分析、模糊决策、模糊评判、系统理论、信息检索、医学、生物学等各个方面。在气象、结构力学、控制、心理 学等方面已有具体的研究成果。然而模糊数学最重要的应用领域是计算机职能,不少人认为它与新一代计算机的研制有密切的联系。
目前,世界上 发达国家正积极研究、试制具有智能化的模糊计算机,1986年日本山川烈博士首次试制成功模糊推理机,它的推理速度是1000万次/秒。1988年,我国 汪培庄教授指导的几位博士也研制成功一台模糊推理机——分立元件样机,它的推理速度为1500万次/秒。这表明我国在突破模糊信息处理难关方面迈出了重要 的一步。
模糊数学还远没有成熟,对它也还存在着不同的意见和看法,有待实践去检验。
四:模糊数学的主要应用
1 模糊数学自身的理论研究进展迅速。我国模糊数学自身的理论研究仍占模糊数学及其应用学科的主导地位,所取得的研究成果在《模糊数学》、《模糊系统与数学》 等数十种学术期刊和全国高校学报中经常可见,模糊聚类分析理论、模糊神经网络理论和各种新的模糊定理及算法不断取得进展。
2.模糊数学目前在自动控制技术领域仍然得到最广泛的应用,所涉及的技术复杂繁多,从微观到宏观、从地下到太空无所不有,在机器人实时控制、电磁元件自适应控制、各种物理及力学参数反馈控制、逻辑控制等高新技术中均成功地应用了模糊数学理论和方法。
3.模糊数学在计算机仿真技术、多媒体辨识等领域的应用取得突破性进展,如图像和文字的自动辨识、自动学习机、人工智能、音频信号辨识与处理等领域均借助了模糊数学的基本原理和方法。
4. 模糊聚类分析理论和模糊综合评判原理等更多地被应用于经济管理、环境科学、安全与劳动保护等领域,如房地价格、期货交易、股市情报、资产评估、工程质量分 析、产品质量管理、可行性研究、人机工程设计、环境质量评价、资源综合评价、各种危险性预测与评价、灾害探测等均成功地应用了模糊数学的原理和方法。
5.地 矿、冶金、建筑等传统行业在处理复杂不确定性问题中也成功地应用了模糊数学的原理和方法,从而使过去凭经验和类比法等处理工程问题的传统做法转向数学化、 科学化,如矿床预测、矿体边界确定、油水气层的识别、采矿方法设计参数选择、冶炼工艺自动控制与优化、建筑物结构设计等都有应用模糊数学的成功实践。
6. 我国医药、生物、农业、文化教育、体育等过去看似与数学无缘的学科也开始应用模糊数学的原理和方法,如计算机模糊综合诊断、传染病控制与评估、人体心理及 生理特点分析、家禽孵养、农作物品种选择与种植、教学质量评估、语言词义查找、翻译辨识等均有一些应用模糊数学的实践,并取得很好效果。
模糊数学目前在自动控制技术领域仍然得到最广泛的应用,所涉及的技术复杂繁多,从微观到宏观、从地下到太空无所不有,在机器人实时控制、电磁元件自适应控制、各种物理及力学参数反馈控制、逻辑控制等高新技术中均成功地应用了模糊数学理论和方法。
模糊数学在计算机仿真技术、多媒体辨识等领域的应用取得突破性进展,如图像和文字的自动辨识、自动学习机、人工智能、音频信号辨识与处理等领域均借助了模糊数学的基本原理和方法。
模 糊聚类分析理论和模糊综合评判原理等更多地被应用于经济管理、环境科学、安全与劳动保护等领域,如房地价格、期货交易、股市情报、资产评估、工程质量分 析、产品质量管理、可行性研究、人机工程设计、环境质量评价、资源综合评价、各种危险性预测与评价、灾害探测等均成功地应用了模糊数学的原理和方法。
地 矿、冶金、建筑等传统行业在处理复杂不确定性问题中也成功地应用了模糊数学的原理和方法,从而使过去凭经验和类比法等处理工程问题的传统做法转向数学化、 科学化,如矿床预测、矿体边界确定、油水气层的识别、采矿方法设计参数选择、冶炼工艺自动控制与优化、建筑物结构设计等都有应用模糊数学的成功实践。
我 国医药、生物、农业、文化教育、体育等过去看似与数学无缘的学科也开始应用模糊数学的原理和方法,如计算机模糊综合诊断、传染病控制与评估、人体心理及生 理特点分析、家禽孵养、农作物品种选择与种植、教学质量评估、语言词义查找、翻译辨识等均有一些应用模糊数学的实践,并取得很好效果。
五:
但是一些概率论 学者认为模糊数学不过是概率论的一个应用而已。一些搞理论数学的人说这不是数学。搞应用的人则说道理说的很好,但真正的实际效果没有。然而,国际著名的应 用数学家考夫曼(A.Kauffman)教授在访华时说:“他们的攻击是毫无道理的,不必管人家说什么,我们努力去做就是。”
模糊数学是一门崭新的数学学科,它的产生不仅拓广了经典数学的基础,而且是使计算机科 学向人们的自然机理方面发展的重大突破。它在科学技术、经济发展和社会学等问题的广泛应用领域中显示了巨大的力量。它虽然只有二十多年的历史,但已被国内 外数学界以及信息、系统、计算机和自动控制科学、人员的普遍关注,它是正在迅速发展中的有着广阔应用前景的一门崭新学科。

2008年6月17日星期二

CamShift算法,OpenCV实现1--Back Projection

CamShift算法,即"Continuously Apative Mean-Shift"算法,是一种运动跟踪算法。它主要通过视频图像中运动物体的颜色信息来达到跟踪的目的。我把这个算法分解成三个部分,便于理解:
1) Back Projection计算
2) Mean Shift算法
3) CamShift算法
在这里主要讨论Back Projection,在随后的文章中继续讨论后面两个算法。

Back Projection
计算Back Projection的步骤是这样的:
1. 计算被跟踪目标的色彩直方图。在各种色彩空间中,只有HSI空间(或与HSI类似的色彩空间)中的H分量可以表示颜色信息。所以在具体的计算过程中,首先将其他的色彩空间的值转化到HSI空间,然后会其中的H分量做1D直方图计算。
2. 根据获得的色彩直方图将原始图像转化成色彩概率分布图像,这个过程就被称作"Back Projection"。
在OpenCV中的直方图函数中,包含Back Projection的函数,函数原型是:
void cvCalcBackProject(IplImage** img, CvArr** backproject, const CvHistogram* hist);
传递给这个函数的参数有三个:
1. IplImage** img:存放原始图像,输入。
2. CvArr** backproject:存放Back Projection结果,输出。
3. CvHistogram* hist:存放直方图,输入

下面就给出计算Back Projection的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 );
int hist_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.计算Back Projection:
IplImage* rawImage;
//----------------------------------------------
//get from video frame,unsigned byte,one channel
//----------------------------------------------
IplImage* result=cvCreateImage(cvGetSize(rawImage),IPL_DEPTH_8U,1);
cvCalcBackProject(&rawImage,result,hist);
5.结果:result即为我们需要的.

算法分析
用 在cvCalcBackProject处理中的模板是目标图像色调(HUE)的直方图,而直方图可以看作是一种概率分布图。在处理前,目标图像中的每一个 象素的值描述的在这一点的颜色信息,而处理后,图像中每一个象素的值就变成了这个颜色信息出现在此处的可能性的一种离散化的度量,出现的可能性大,象素的 值就大,反之则小。这样就为后面的匹配和跟踪提供了线索。

2008年6月4日星期三

Kalman filter

1 什么是卡尔曼滤波器
What is the Kalman Filter?

在学习卡尔曼滤波器之前,首先看看为什么叫卡尔曼。跟其他著名的理论(例如傅立叶变换,泰勒级数等等)一样,卡尔曼也是一个人的名字,而跟他们不同的是,他是个现代人!

卡尔曼全名Rudolf Emil Kalman,匈牙利数学家,1930年出生于匈牙利首都布达佩斯。19531954年于麻省理工学院分别获得电机工程学士及硕士学位。1957年于哥伦比亚大学获得博士学位。我们现在要学习的卡尔曼滤波器,正是源于他的博士论文和1960年发表的论文《A New Approach to Linear Filtering and Prediction Problems》(线性滤波与预测问题的新方法)。如果对这编论文有兴趣,可以到这里的地址下载: http://www.cs.unc.edu/~welch/media/pdf/Kalman1960.pdf

简单来说,卡尔曼滤波器是一个“optimal recursive data processing algorithm(最优化自回归数据处理算法)。对于解决很大部分的问题,他是最优,效率最高甚至是最有用的。他的广泛应用已经超过30年,包括机器人导航,控制,传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理,例如头脸识别,图像分割,图像边缘检测等等。

2
.卡尔曼滤波器的介绍
Introduction to the Kalman Filter

为了可以更加容易的理解卡尔曼滤波器,这里会应用形象的描述方法来讲解,而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是,他的5条公式是其核心内容。结合现代的计算机,其实卡尔曼的程序相当的简单,只要你理解了他的那5条公式。

在介绍他的5条公式之前,先让我们来根据下面的例子一步一步的探索。

假设我们要研究的对象是一个房间的温度。根据你的经验判断,这个房间的温度是恒定的,也就是下一分钟的温度等于现在这一分钟的温度(假设我们用一分钟来做时间单位)。假设你对你的经验不是100%的相信,可能会有上下偏差几度。我们把这些偏差看成是高斯白噪声(White Gaussian Noise),也就是这些偏差跟前后时间是没有关系的而且符合高斯分配(Gaussian Distribution)。另外,我们在房间里放一个温度计,但是这个温度计也不准确的,测量值会比实际值偏差。我们也把这些偏差看成是高斯白噪声。

好了,现在对于某一分钟我们有两个有关于该房间的温度值:你根据经验的预测值(系统的预测值)和温度计的值(测量值)。下面我们要用这两个值结合他们各自的噪声来估算出房间的实际温度值。

假如我们要估算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,就是卡尔曼增益(Kalman Gain)。他可以随不同的时刻而改变他自己的值,是不是很神奇!

下面就要言归正传,讨论真正工程系统上的卡尔曼。

3
卡尔曼滤波器算法
The Kalman Filter Algorithm

在这一部分,我们就来描述源于Dr Kalman 的卡尔曼滤波器。下面的描述,会涉及一些基本的概念知识,包括概率(Probability),随即变量(Random Variable),高斯或正态分配(Gaussian Distribution)还有State-space Model等等。但对于卡尔曼滤波器的详细证明,这里不能一一描述。

首先,我们先要引入一个离散控制过程的系统。该系统可用一个线性随机微分方程(Linear Stochastic Difference equation)来描述:
X(k)=A X(k-1)+B U(k)+W(k)
再加上系统的测量值:
Z(k)=H X(k)+V(k)
上两式子中,X(k)k时刻的系统状态,U(k)k时刻对系统的控制量。AB是系统参数,对于多模型系统,他们为矩阵。Z(k)k时刻的测量值,H是测量系统的参数,对于多测量系统H为矩阵。W(k)V(k)分别表示过程和测量的噪声。他们被假设成高斯白噪声(White Gaussian Noise),他们的covariance 分别是QR(这里我们假设他们不随系统状态变化而变化)。

对于满足上面的条件(线性随机微分系统,过程和测量都是高斯白噪声),卡尔曼滤波器是最优的信息处理器。下面我们来用他们结合他们的covariances 来估算系统的最优化输出(类似上一节那个温度的例子)。

首先我们要利用系统的过程模型,来预测下一状态的系统。假设现在的系统状态是k,根据系统的模型,可以基于系统的上一状态而预测出现在状态:
X(k|k-1)=A X(k-1|k-1)+B U(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)=A P(k-1|k-1) A’+Q ……… (2)
(2)中,P(k|k-1)X(k|k-1)对应的covarianceP(k-1|k-1)X(k-1|k-1)对应的covarianceA’表示A的转置矩阵,Q系统过程的covariance。式子12就是卡尔曼滤波器5个公式当中的前两个,也就是对系统的预测。

现在我们有了现在状态的预测结果,然后我们再收集现在状态的测量值。结合预测值和测量值,我们可以得到现在状态(k)的最优化估算值X(k|k)
X(k|k)= X(k|k-1)+Kg(k) (Z(k)-H X(k|k-1)) ……… (3)
其中Kg为卡尔曼增益(Kalman Gain)
Kg(k)= P(k|k-1) H’ / (H P(k|k-1) H’ + R) ……… (4)

到现在为止,我们已经得到了k状态下最优的估算值X(k|k)。但是为了要另卡尔曼滤波器不断的运行下去直到系统过程结束,我们还要更新k状态下X(k|k)covariance
P(k|k)=
I-Kg(k) HP(k|k-1) ……… (5)
其中I 1的矩阵,对于单模型单测量,I=1。当系统进入k+1状态时,P(k|k)就是式子(2)P(k-1|k-1)。这样,算法就可以自回归的运算下去。

卡尔曼滤波器的原理基本描述了,式子12345就是他的5 个基本公式。根据这5个公式,可以很容易的实现计算机的程序。

下面,我会用程序举一个实际运行的例子。。。
4
简单例子
A Simple Example

这里我们结合第二第三节,举一个非常简单的例子来说明卡尔曼滤波器的工作过程。所举的例子是进一步描述第二节的例子,而且还会配以程序模拟结果。

根据第二节的描述,把房间看成一个系统,然后对这个系统建模。当然,我们见的模型不需要非常地精确。我们所知道的这个房间的温度是跟前一时刻的温度相同的,所以A=1。没有控制量,所以U(k)=0。因此得出:
X(k|k-1)=X(k-1|k-1) ……….. (6)
式子(2)可以改成:
P(k|k-1)=P(k-1|k-1) +Q ……… (7)

因为测量的值是温度计的,跟温度直接对应,所以H=1。式子345可以改成以下:
X(k|k)= X(k|k-1)+Kg(k) (Z(k)-X(k|k-1)) ……… (8)
Kg(k)= P(k|k-1) / (P(k|k-1) + R) ……… (9)
P(k|k)=
1-Kg(k)P(k|k-1) ……… (10)

现在我们模拟一组测量值作为输入。假设房间的真实温度为25度,我模拟了200个测量值,这些测量值的平均值为25度,但是加入了标准偏差为几度的高斯白噪声(在图中为蓝线)。

为了令卡尔曼滤波器开始工作,我们需要告诉卡尔曼两个零时刻的初始值,是X(0|0)P(0|0)。他们的值不用太在意,随便给一个就可以了,因为随着卡尔曼的工作,X会逐渐的收敛。但是对于P,一般不要取0,因为这样可能会令卡尔曼完全相信你给定的X(0|0)系统最优的,从而使算法不能收敛。我选了X(0|0)=1度,P(0|0)=10

系统的真实温度为25度,图中用黑线表示。图中红线是卡尔曼滤波器输出的最优化结果(该结果在算法中设置了Q=1e-6R=1e-1)。