跟踪与运动
TODO
第10章
跟踪与运动
跟踪基础
当我们在处理一段视频而非某张静止的图像时，通常视频中有一个或几个特定物体
是我们想在整个视野范围内关注的。在前一章中，我们知道了如何在逐帧（frame-
by-frame)分离出特定的形状，例如人或者车；现在，我们尝试理解这个物体的运
动。理解物体的运动则主要包含两个部分：识别（identification)和建模。
识别指在视频流后续的帧中找出之前某帧中的感兴趣物体。前面章节中讲到的矩和
颜色直方图可以帮助识别我们关注的物体。一个相关的问题就是跟踪不明物体。当
我们需要根据运动来确定什么是感兴趣的，或者当物体的运动使得其成为感兴趣物
体时，跟踪不明物体就很重要了。经典的跟踪不明物体的方法是跟踪视觉上重要的
关键点，并不是整个物体。OpenCV 提供了两种方法实现跟踪关键点：Lucas-
Kanade?[Lucas81]和 Horn-Schunk [Horn81]方法。这两种方法分别代表了通常提到
的稀疏（sparse)和稠密（dense)光流。
上述方法实际上只能给出物体实际位置的初步计算。理解物体运动的第二个部分，
也就是建模，则可以帮助我们解决这个问题。为估计由这种粗略测量得出的物体的
轨迹，人们设计了许多有力的数学方法。这些方法适用于二维和三维物体或物体位
置的模型。
①很奇怪，OpenCV 中实现的金字塔框架下的Lucas-Kanade 光流方法是根据一篇未正式
发表的论文实现的，这篇论文是Bouguet的［Bouguet04]:
寻找角点
有多种局部特征可以用来进行跟踪。花一点时间去思考究竟角点（corner)用什么样
的特征来描述是值得的。很明显，如果从一面很大的空墙上选择一个点，那么将很
难从视频的下一帧中再找出这一点。如果墙上的所有点都是一样的或者是相似的，
我们就不会有太好的运气能在随后的视频帧中跟踪到这个点了。相反，如果选择一
个独一无二的点，那么再找到这点的几率就非常大。实际上，被选择的点或特征应
该是独一无二的，或者至少接近独一无二，并且在与另一张图像的其他点可以进行
参数化的比较。如图10-1所示。
【316~317】
图10-1:圆圈圈出的特征点是易于跟踪的点，而矩形圈出的特征－
甚至是明
显的边缘－却不然
回到那面大而空的墙，我们想要找本身具有明显变化的点－例如导数值比较明显
的地方。虽然仅此是不够的，但这是一个开始。一个导数值比较明显的点可能在某
种类型的缘上，但却和此边缘上的其他点看起来一样（参考图 10-8和“Lucas-
Kanade方法”一节中讨论的孔径问题）。
如果一个点在两个正交的方向上都有明显的导数，则我们认为此点更倾向是独一无
二的，所以，许多可跟踪的特征点都称为角点。从直观上讲，角点（而非边缘）是一
类含有足够信息且能从当前帧和下一帧中都能提取出来的点。
最普遍使用的角点定义是由Harris[Harris88]提出的。定义的基础是图像灰度强度
的二阶导数（x,?2y,dxdy)矩阵。考虑到图像所有的像素点，我们可以想像图像的二
阶导数即形成一幅新的“二阶导数图像”，或融合在一起形成一幅新的 Hessian 图
像。这个术语源自于一个点的如下定义的二维Hessian矩阵：
跟踪与运动
H(p)=
【317】
对于 Harris 角点，我们使用每点周围小窗口的二阶导数图像的自相关矩阵。这个
自相关矩阵的定义如下：
2(x+iy+j)
（x+i,y+j),(x+y+17
M(x,y)=
3
w.1(x+i,y+j)I,(x+i,y+j)
（x+3y+)
－KSI,jsK
这里w,是可以归一化的权重比例，但是通常被用作产生圆形窗口或高斯权重。
Harris 定义的角点位于图像二阶导数的自相关矩阵有两个最大特征值的地方，这在
本质上表示以此点为中心周围存在至少两个不同方向的纹理（或者边缘），正如实际
的角点是由至少两个边缘相交于一点而产生。之所以采用二阶导数是由于它对均匀
梯度不产生响应?.角点的这个定义还有另一个优点。被跟踪的物体在移动过程中
也可能会旋转，找到同时对移动和旋转不变的量是很重要的。只考虑自相关矩阵的
特征值可达到这个目的。这两个最大特征值不仅可以判定一个点是否一个好的可跟
踪的特征点，同时也提供了对这个点进行识别的一个标识。
Harris 最原始的定义是将矩阵H(p)的行列式值与H(p)的迹（带权重系数）相减，再
将差值同预先给定的阈值进行比较。后来 Shi 和Tomasi[Shi94]发现，若两个特征
值中较小的一个大于最小阈值，则会得到强角点。Shi 和 Tomasi 的方法比较充
分，并且在很多情况可以得到比使用Harris方法更好的结果。
函数cvGoodFeaturesToTrack()采用 Shi和 Tomasi提出的方法，先计算二阶导数
（利用Sobel算子），再计算特征值，它返回满足易于跟踪的定义的一系列点。
void cvGoodFeaturesToTrack(
const CvArr*
image,
CvArr*
eigImage,
CvArr*
tempImage,
CvPoint2D32f*
corners,
int*
corner_count,
double
quality_level,
double
min_distance,
①
梯度是一阶导数而来。若一阶导数是均匀的（常数），则二阶导数为0.
st
CvArr*
mask
＝NULL,
block_size
＝3,
use_harris
＝0,
ouble
k
＝0.4
【318】
数中，
输入图像 image须为8位或32位（也就是，IPL_DEPTH_8U 或者
32F)单通道图像。第二和第三个参数是大小与输人图像相同的32位单
参数
tempImage 和 eigImage 在计算过程中被当作临时变量使用，计
量eigImage中的内容是有效的。特别地，每个元素包含了输入图像中
最小特征值。
corners 是函数的输出，为检测到的32位（CvPoint2D32f)
组，在调用函数cvGoodFeaturesTorrack()之前需要给这个数组分配内
很明显，只能为此数组分配有限的内存。corner_count表示可以返回的
数目。函数调用结束后，corner_count 输出实际检测到的角点数目。
level
表示一点被认为是角点的可接受的最小特征值。实际用于过滤角
最小特征值是
quality_level 与图像中最大特征值的乘积。所以
1eve1的值不应超过1(常用的值为0.10或者0.01).检测完所有的角点
要进一步剔除掉一些距离较近的角点。min_distance 保证返回的角点之间
5不小于min_distance个像素。
是可选参数，是一幅像素值为布尔类型的图像，用于指定输人图像中参与角
的像素点。若mask的值设为NOLL,则选择整个图像。block_size是计算
自相关矩阵时指定点的领域，采用小窗口计算的结果比单点（也就是
size为1)计算的结果要好。若use_harris的值为非0,则函数使用
的角点定义；若为0,则使用shi-Tomasi的定义。当 use_harris为k且
则k为用于设置Hessian 自相关矩阵即对 Hessian行列式的相对权重的权重
cvGoodreaturesTorrack()的输出结果为角点的位置数组，我们往往会希望
第一幅相似的图像中寻找这些角点。就本章而言，我们关注的是在后续的视频帧
我到这些特征点，但是这个函数也有其他的应用，例如，对从微小不同的视角拍
多幅图像进行匹配。这个问题将在后面讨论立体视觉的小节中遇到。
花
亚像素级角点
如果我们进行图像处理的目的不是提取用于识别的特征点而是进行几何测量，则通
常需要更高的精度，而函数cVGoodFeaturesToTrack()只能提供简单的像素的坐
标值，就是说，我们有时会需要实数坐标值而不是整数坐标值，例如（8.25,
117.16).
您也许需要确定图像中一个尖锐的峰值点的位置，但是峰值的位置一般都不会恰好
位于一个像素的正中心，使得确定其位置比较困难。要解决这个问题，需要先用一
条曲线（例如一条抛物线）拟合图像的值，再用数学的方法计算位于像素之间的峰值
点的位置。亚像素检测方法就是一些有关峰值点位置的计算技巧（了解更多此方面
知识以及最新的方法，请参考Lucchese[Lucchese02]和Chen[Chen05]).图像测
量常用的领域为三维重建、摄像机标定、图像拼接以及在卫星图像中查找特定信
号，如一栋建筑的精确位置。
【319~320】
亚像素级角点的位置在摄像机标定、跟踪并重建摄像机的轨迹或者重建被跟踪目标
的三维结构时是一个基本的测量值。我们已经知道如何对整数网格像素来检测角点
并得到角点的位置，下面所要讨论的是如何将求得的角点位置精确到亚像素级精
度。一个向量和与其正交的向量的点积为0,角点则满足这种情况，如图10-2
所示。
邻域
边上的梯度
VI(P)
0=(d)I
（a):p点在一个区域内部
（b):p点在区域边缘
在这两种情况下，点积都为零
＜VI(p),q-p>=0
图10-2:计算亚像素级精度的角点：（a)点p附近的图像是均匀的，其梯度为
0;(b)边缘的梯度与沿边缘方向的q-p向量正交。在图中两种情况下，p点的
梯度与q-p向量的点积均为0(参考正文）
图10-2假设起始角点q在实际亚像素级角点的附近。检测所有的q-p向量。若点p
位于一个均匀的区域，则点p处的梯度为0.若q-p 向量的方向与边缘的方向一
致，则此边缘上p点处的梯度与q-p向量正交，在这两种情况下，p点处的梯度与
q-p向量的点积为0.我们可以在p点周围找到很多组梯度以及相关的向量q-p,令
其点集为0,然后可以通过求解方程组，方程组的解即为角点q的亚像素级精度的
位置，也就是精确的角点位置。
【320】
函数cvFindCornerSubPix()用于发现亚像素精度的角点位置。
void cvFindCornerSubPix(
const CvArr*
image,
CvPoint2D32f*
corners,
int
count,
CvSize
win,
CvSize
zero_zone,
CvTermCriteria
criteria
）；
输入图像 image是8位单通道的灰度图像。corners 为整数值的像素位置，例如
是由cVGoodFeatures()得到的位置坐标。corners
设定了角点的初始位置。
count为需要计算的角点数目。
实际计算亚像素级的角点位置时，解的是一个点积的表达式为0的方程组（参考图
10-2),其中每一个方程都是由q邻域的一个点产生。win 指定了等式产生的窗口
的尺寸。搜索窗口的中心是整数坐标值的角点并从中心点在每个方向上扩展 win
指定的像素（也就是说，如果win.width=4,则实际的窗口大小为4+1+4=9个像
素宽）。这些等式构成一个可用自相关矩阵的逆（非前文讨论 Harris 角点时涉及的
自相关矩阵）来求解的线形方程组。实际上，由于非常接近p的像素产生了很小的
特征值，所以这个自相关矩阵并不总是可逆的。为了解决这个问题，一般可以简单
地剔除离p点非常近的像素。输入参数zero_zone定义了一个禁区（与win相似，
但通常比 win小），这个区域在方程组以及自相关矩阵中不被考虑。如果不需要这
样一个禁区，则zero_zone应设置为cvSize(-1,-1).
当找到一个q的新位置时，算法会以这个新的角点作为初始点进行迭代直到满足用
户定义的迭代终止条件。迭代过程的终止条件可以是最大迭代次数
CV_TERMCRIT_ITER 类型，或者是设定的精度 CV_TERMCRIT_ITER 类型（或者是两
者的组合）。终止条件的设置在极大程度上影响最终得到的亚像素值的精度。例
如，如果指定0.10,则求得的亚像素级精度为像素的十分之一。
不变特征
继 Harris 提出角点以及后来 Shi和 Tomasi提出角点后，许多其他类型的角点和相
关局部特征点也被提出来。SIFT(“scale-invariant feature transform”）是其中一种
广泛应用的类型［Lowe04].SIFT 特征正如其名称一样是缩放不变的。SIFT 在一点
处检测主要梯度方向，根据这个方向记录局部梯度直方图结果，所以 SIFT 也是旋
转不变的。正是由于这些特点，SIFT 特征在小的仿射变换中有相对不错的表现。
虽然OpenCV库中还没有实现SIFT算法（参考第14章），但可以用OpenCV的基本
函数实现它。我们将不在这个问题上花更多的时间，但需要记住的是，用前面已经
介绍过的 OpenCV 函数可以实现计算机视觉领域的文献中提出的绝大部分特征。
【321~322】
光流
如前所述，我们经常不知道任何关于视频内容的先验知识，但是需要估计两帧之间
（或一个帧序列）的运动。运动本身就说明了存在运动的感兴趣的目标。在图10-3
中解释了光流。
a
光流
2
23
I(t+1)
（a)/(t),点集
：｛pi}
（b)速度向量｛3}
图像跟踪
（c)图像序列（单个摄像机）
（d)跟踪结果
图10-3:光流。目标特征点（a)被连续跟踪，它们的运动被转换成速度矢量
（b);下面的两幅图分别表示ー幅走廊的图像（c)和摄像机沿着走廊运动时的光
流矢量（d)(原图由 Jean-Yves Bouguet提供）
可以将图像中的每个像素与速度关联，或者等价地，与表示像素在连续两帧之间的
位移关联。这样得到的是稠密光流（dense optical flow),即将图像中的每个像素都
数
与速度关联。Horm-Schunck 方法［Horn81]计算的就是稠密光流的速度场。OpenCV
中实现了一种比较简单直接的方法，即对前后连续两帧的一个像素的邻域进行匹
配。这种方法被称为块匹配（block matching).这两种方法都会在“稠密跟踪方法”
一节中讨论。
【322~323】
实际上计算稠密光流并不容易。我们先来看看一张白纸的运动。上一帧中的白色的
像素在下一帧中仍然为色，而只有边缘的像素并且是与运动方向垂直的像素才会
产生变化。稠密光流的方法需要使用某种插值方法在比较容易跟踪的像素之间进行
插值以解决那些运动不明确的像素，从中可以清楚地看到稠密光流相当大的计
算量。
于是我们找到了另一种方法，稀疏光流（sparse optical flow).稀疏光流的计算需要
在被跟踪之前指定一组点。如果这些点具有某种明显的特性，例如前面讲到的“角
点”，那么跟踪就会相对稳定和可靠。OpenCV 可以帮助我们找到最适合跟踪的特
征点。在很多的实际应用中，稀疏跟踪的计算开销比稠密跟踪小得多，以至于后者
只有在理论中研究。
下面的几节介绍了几种不同的跟踪的方法。我们从最流行的稀疏跟踪方法 Lucas-
Kanade(LK)光流开始介绍；这种方法与图像金字塔一起，可以跟踪更快的运动。
接下来介绍两种稠密跟踪方法，Horn-Schunck方法和块匹配方法。
Lucas-Kanade 方法
Lucas-Kanade 算法（LK)[Lucas81]是用于求稠密光流的，最初于1981年提出。由于
算法易于应用在输入图像中的一组点上，其后来成为求稀疏光流的一种重要方法。
LK 算法只需要每个感兴趣点周围小窗口的局部信息，所以它可以应用于稀疏内
容，这与 Horn 和 Schunck 的算法全局性是不同的（后面会详细讲述）。但是，使用
小窗口的LK算法存在不足之处，较大的运动会将点移出这个小窗口，从而造成算
法无法再找到这些点。金字塔的LK算法可以解决这个问题，即从图像金字塔的最
高层（细节最少）开始向金字塔的低层（丰富的细节）进行跟踪。跟踪图像金字塔允许
小窗口捕获较大的运动。
这个算法是一个重要而有效的方法，我们将会具体介绍其中的数学问题。不想阅读
①
Black 和 Anadan 提出的计算稠密光流的方法［Black93;Black96]被广泛用在电影制作
中，因为为了追求视觉的质量，电影工作室愿意牺牲时间开销来获得详细的光流信
息。OpenCV未来的版本中会将这些方法包含进去（参考第14章）。
这些细节的读者可以直接阅读函数描述和代码部分，但是建议读者至少浏览一遍其
中描述Lucas-Kanade 光流假设条件的文字和插图部分，这样在失去头绪的时候知
道应该如何去解决。
Lucas-Kanade 算法原理
LK算法基于以下三个假设。
（1)亮度恒定。图像场景中目标的像素在帧间运动时外观上保持不变。对于灰度图
像（LK算法也可用于彩色图像），需要假设像素被逐帧跟踪时其亮度不发生变化。
（2)时间连续或者运动是“小运动”。图像的运动随时间的变化比较缓慢。实际应
用中指的是时间变化相对图像中运动的比例要足够小，这样目标在帧间的运动
就比较小。
（3)空间一致。一个场景中同一表面上邻近的点具有相似的运动，在图像平面上的
投影也在邻近区域。
图10-4 说明了这三点假设如何构建一个有效的跟踪算法。第一点假设：亮度恒
定，指对被跟踪部分（patch)像素不随时间变化的要求：
（p+1*(p+1)x)I=(r(x)x)I=(r*x)f
假设
亮度恒定
（a)
1(x+u,y+v,t+1)=1(x,y,t)
时间连续
空间一致
表面
（q)
（c)
图像平面
图10-4:Lucas-Kanade 光流假设：（a)场景中物体被跟踪的部分的亮度不变；
（b)运动相对于帧率是缓慢的；（c)[相邻的点保持相邻（图由 Michael Black
［Black82]提供）］
【324】
很简单，公式的意思是被跟踪像素的灰度不随时间变化而变化：
af(x)=0
10
第二点假设：时间连续，指相邻帧之间的运动较小。换句话，可以将运动的变化看
成是亮度对时间的导数（即断言一个序列中相邻帧间的变化是小微分的）。我们先来
看一下一维空间的例子，以理解第二点假设所蕴含的意思。
由亮度恒定的等式开始，考虑隐含的x为t的函数这点，将亮度的定义f(x,t)用
I(x(t),t)替换，再应用偏微分的链式规则（chain rule),即：
o-
其中，I,是图像的偏导数，I,是图像随时间的导数，v是要求的速度。简单一维空
间中的光流速度的等式为：
＝－1/1
下面我们对一维跟踪问题进行一些直观的解释。图 10-5 表示的是左边值大右边值
小、沿x轴向右运动的一个“边缘”。我们要求出图的上半部分标出的边缘运动的
速度V.在图的下半部分，对这个速度的测量正是“rise over run”，即 rise 是随
时间变化的而run是斜率（空间的导数）。负号纠正了x的斜率。
一维光流
1(x,t)
［（x,t+1)
I(x,t)
1(x,t+1)
时域导数
空域导数
假定：
?亮度不变
／e=
＊小运动
图10-5:一维空间中的Lucas-Kanade光流。通过计算亮度对时域导数与亮度
对空域导数的比值，可以估计出运动边缘（图的上半部分）的速度
图10-5 揭示了光流公式的另一个问题：我们的假设可能不是十分正确，即图像的
亮度实际上不是恒定的，时间步长（由摄像机设定）也不是如我们期望的相对于运动
足够短。所以我们求解的速度并不准确。但是，如果求解到的速度与实际的速度
“足够接近”，就可以用迭代的方法来解决这个问题，如图 10-6所示，将第一次
估计的速度（不准确的）作为初始值进行下一次迭代并重复这个过程。根据沿x运动
的像素不变的亮度恒定假设，在迭代过程中，可以保持使用由第一帧计算得到的x
的空间导数。这种重复使用已经计算出的空间导数极大地节省了计算开销。时间导
数仍需要在每次迭代和每一帧中重新计算，但如果初始值与实际值足够接近，迭代
过程在5次迭代之内就会收敛到接近的准确值，这就是所谓的牛顿法。如果初始估
计值与实际值相差比较大，则牛顿法实际上是发散的。
【325~326】
代可以优化速度向量
I(x,t)
I(x,t+1)
第2次迭代得到的时域导数
x
J
空域导数也采用
Fpmerion 4/4
约5次迭代后收敛
图10-6:通过迭代来优化光流的求解（牛顿法）。使用同样的两个图像和同样的
空域导数，重新求解时域导数，通常可通过几次迭代收敛到稳定解
前面是LK算法在一维空间的用法，下面我们将其扩展应用到二维空间的图像上。
乍看起来好像比较简单，只需要加入y坐标即可。改变一下符号，速度的y分量为
v,x分量为u,则得到：
Iu+I,v+I,=0
不幸的是，这个等式对任意一个像素都含有两个未知量。这就是说对于单个像素，
等式的约束条件过少，不能得到此点二维运动的定解。我们只能求得与光流方程描
述的方向垂直的运动分量。图10-7给出了数学和几何的细节。
垂直光流由孔径问题产生，即用小孔或小窗口去测量运动。这种情况下，我们通常
只能观测到边缘而观测不到角点，而只依靠边缘是不足以判断整个物体是如何运动
（也就是朝哪个方向运动）的，如图10-8所示。
6九
从一维到二维光流
公式
但是，只能求得一条直线，而不是一个点
Ixu+Iyv+I4=0
Ixu+Iyv=-14
I-=n IΔ
1,u+1,+1=0
＝1A=
＂Normal flow”
图10-7:单个像素的二维光流：单个像素的光流是无法求解的，最多只能求出
光流方向，因为光流方向与光流方程描述的线垂直（图由Michael Black提供）
（a)
孔经问题
（b)
图10-8:孔径问题：从aperture window(a)我们可以观测到边缘向右运动，但
是无法观察到边缘也在向下运动（b)
那么，如何解决单个像素不能求解整个运动的问题呢？这时需要利用光流的最后一
点假设。若一个局部区域的像素运动是一致的，则可以建立邻域像素的系统方程来
求解中心像素的运动。例如，如果用当前像素 5x5①邻域的像素的亮度值（彩色像素
的光流只需增加两倍）来计算此像素的运动，则可以建立如下的25个方程。
【326~327】
①
窗口可以是 3x3、7x7 或者任何被指定的值。若窗口太大则会由于违背运动一致的假
设而不能进行较好的跟踪。若窗口太小，则又会产生孔径问题。
（
1,(p1)
I,(P2) I,(P2)
I,(P2s)
1,(P25)
25x2
现在我们得到一个约束条件过多的系统方程，若在5x5的窗口中包含两条或院上
边缘则可以解此系统方程。为了解这个系统方程，需要建立一个该方程的最小事
方，通过下面方程来求解最小化的｜Ad-bP:
／＝（
由这个关系式可以得到u和v运动分量。这个关系更详细的表述如下：
S11, ΣI1,u
ΣIIy
Σ1,
当（AA)可逆时，方程的解如下：
＝（AA)-1ATb
【328】
当（AA)满秩（秩为2)也即（AA)有两个较大特征向量时，（AA)可逆。图像中纹理至
少有两个方向的区域，这个条件可以满足。这种情况下，跟踪窗口的中心在图像的
角点区域时，（AA)的特性最好。由这一点我们可以联系到前面讨论的 Harris 角点
检测器。事实上，正因为（AA)在角点出有两个大的特征向量，所以这些角点是
“可用于跟踪的良好特征点”（参考前面对 cvGoodFeaturesToTrack()的介绍），
后文即将介绍的cvCalcopticalFlowLK()实现了这个计算。
理解了小而连贯运动的假设的读者现在会感到不解：对于大多数30Hz的接像机。
大而不连贯的运动是普遍存在的情况。而 Lucas-Kanade 光流正因为这个原因在实
际中的跟踪效果并不是很好：我们需要一个大的窗口来捕获大的运动，面大窗口往
往会违背运动连贯的假设！图像金字塔可以解决这个问题，即最初有铰大的空间尺
度上进行跟踪，再通过对图像金字塔向下直到图像像素的处理来修正初始运动速度
的假定。
所以，建议的跟踪方法是：在图像金字塔的最高层计算光流，用得到的运动基计结
果作为下一层金字塔的起始点，重复这个过程直到到达金字塔的最底层。这样就将
不满足运动假设的可能性降到最小从而实现对更快和更长的运动的跟踪。这个更精
致的算法叫金字塔Lucas-Kanade 光流，其原理在图10-9中描述。实现金字塔
Lucas-Kanade 光流的函数是 cvCalcOpticalFlowPyrLK(),后面将介绍这个函数。
从粗到精的光流估计
run iterative L-K
warp and upsample
run iterative L-K<
image1
高斯金字塔图像
Ic-1
高斯金字塔图像
图 10-9:金字塔Lucas-Kanade光流。为了减轻由小而连贯的运动假设引起的
问题，首先在金字塔顶层计算光流；在上一次估计到的运动将作为下一层的起
始点，进行进一步估计
Lucas-Kanade 代码
下面的函数实现了非金字塔的Lucas-Kanade稠密光流算法：
void cvCalcopticalFlowLK(.
const CvArr*
imgA,
const CvArr*
imgB,
CvSize
winSize,
CvArr*
velx,
CvArr*
vely
）；
OpenCV 这个函数的输出只记录了可以计算最小误差的像素。对于最小误差不能被
可靠地计算出的像素，相关的速度被设置为0.多数情况下我们不会使用这个函
数，而用下面介绍的基于图像金字塔的方法。
金字塔Lucas-Kanande 代码
现在来介绍 OpenCV 中在图像金字塔中计算 Lucas-Kanade 光流的算法
cvCalcopticalFlowPyrLK().我们将会看到，这个计算光流的函数使用了“易于
跟踪的特征点”并返回每个点被跟踪的情况。
【329~330】
void cvCalcopticalFlowPyrLK(
const CvArr*
imgA,
const CvArr*
imgB,
CvArr*
pyrA,
CvArr*
pyrB,
CvPoint2D32f*
featuresA,
CvPoint2D32f*
featuresB,
int
count,
CvSize
winSize,
int
level,
char*
status,
float*
track_error,
CvTermCriteria criteria,
int
flags
）；
此函数有很多输入参数，我们现花一点时间来弄清楚它们是什么。一旦掌握了这个
函数，我们就可以到下一个问题，即跟踪哪些点和怎样计算。
cvCalcopticalFlowPyrLk()的前两个参数代表初始图像和最终图像，两幅图像都
是8位的单通道图像。第三、四个参数是申请存放两幅输入图像（pyrA和 pyrB)金
字塔的缓存。这两个缓存的大小至少为（img.width+8)*img.height/3个字节?.
（如果这两个指针设置为空，则当函数被调用时会自动分配合适的内存，使用并释
放之，但这样会降低函数的执行效率。）featuresA 数组存放的是用于寻找运动的
点，featuresB 与 featuresB 相似，存放 featuresA 中点的新位置；count是
featuresA 中点的数目。win_size定义了计算局部连续运动的窗口尺寸。参数
level 用于设置构建的图像金字塔的栈的层数。若 level 设置为0则不使用金字
塔。数组 status的长度是 count;函数调用结束时，status 中的每个元素被置
1(对应点在第二幅图像中被发现）或0(对应点在第二幅图像中未被发现）。参数
track_error 为可选参数，它是表示被跟踪点的原始图像小区域与此点在第二幅
图像的小区域间的差的数组。track_error 可以删除那些局部外观小区域随点的
①分配这个大小的缓存是因为演算空间不仅存放图像本身，还要存放整个金字塔。
运动变化剧烈的点。
我们还需要一个迭代终止的条件，即 criteria.CvTermCriteria 是被许多含有
迭代过程的OpenCV算法使用的一个结构：
cvTermCriteria(
int type, // CV_TERMCRIT_ITER, CV_TERMCRIT_EPS, or both
int max_iter,
double epsilon
）；
我们使用cvTermCriteria()函数来生成需要的这个结构。函数的第一个参数为
CV_TERMCRIT_ITER 或者 CV_TERMCRIT_EPS.这个参数分别告诉算法是根据迭代
次数来终止迭代过程，还是根据收敛误差是否达到一个较小值来终止迭代。后面的
两个参数则设置终止算法的一个或两个准则的值。有两个选项的原因使我们可以设
税市桑CERCRIITERICV_TERMCRITEPS,刻兼市
时迭代终止（大多数的实际代码就是这样设置的）。
最后一个参数 flags允许对函数内部 bookkeeping进行一些细微的控制；它可以设
置为下面的任何一个或全部（使用位OR运算）。
CV_LKFLOW_PYR_A_READY 在调用和存储到
pyrA 之前，先计算第一帧
的金字塔。
CV_LKFLOW_PYR_B_READY 在调用和存储到 pyrB之前，先计算第二帧
的金字塔。
CV_LKFLOW_INITIAL_GUESSES 在函数调用之前，数组 B已包含特征点
的初始坐标值。
【331~332】
这些标志特别是在处理视频时会用到。图像金字塔的计算量较大，所以应该尽可能
避免重复计算金字塔。计算得到的图像对的后面一帧被作为下次计算的图像对的初
始帧。如果调用者为函数分配了缓存（而不是让函数在内部分配），则函数返回时，
每幅图像的金字塔将被存储在其中。告诉函数金字塔已经被计算，则函数不会再进
行重复计算。类似地，如果从前一帧中计算得到点的运动，则在计算这些点在下一
帧中的位置就有了一个好的初始估计。
应用的基本过程很简单：输入图像，在 featuresA 中列出需要跟踪的点，然后调
用函数。函数返回后，检查status
数组以确定哪些点被成功跟踪，再检查
featuresB得到这些点新的位置。
现在我们再来讨论前面未讨论的问题：如何确定哪些特征点是易于被跟踪的。前面
介绍过 OpenCV 函数 cvGoodFeaturesToTrack(),这个函数使用 Shi 和 Tomasi
提出的方法可靠地解决了这个问题。在多数情况下，使用 cvGoodFeatures-
ToTrack()和 cvcalcopticalFlowPyrLK()的组合可以获得很好的结果。当然，读
者可以使用自己的准则来确定被跟踪的点。
下面来看一个简单的使用 cvGoodFeaturesToTrack()和 cvCalcopticalFlow-
PyrLK()的组合的例子（例10-1);也可参考图10-10.
Flow vectors
图10-10:金字塔 Lucas-Kanade稀疏光流。中间的图像是左边图像的下一帧
视频；右边的图像表示计算出的“易于跟踪特征点”的运动（右下图为了增加可
视性在黑色背景中标识出光流矢量）
例10-1:金字塔 Lucas-Kanade 光流代码
／／ Pyramid L-K optical flow example
／／
＃include <cv.h>
＃include <cxcore.h>
＃include <highgui. h>
const int MAX_CORNERS =500;
int main(int argc, char** argv) {
／／ Initialize, load two images from the file system, and
／／ allocate the images and other structures we will need for
／／ results.
／／
IplImage* imgA=cvLoadImage ("image0.jpg",CV_LOAD_IMAGE_GRAYSCALE) ;
IplImage* imgB=cvLoadImage ("image1. jpg",CV_LOAD_IMAGE_GRAYSCALE);
CvSize img_sz = cvGetSize( imgA );
int win_size=10;
IplImage* imgC = cvLoadImage(
＂。。／Data/OpticalFlow1. jpg",
CV_LOAD_IMAGE_UNCHANGED
）；
／／ The first thing we need to do is get the features
／／ we want to track.
／／
IplImage* eig_image = cvCreateImage ( img_sz,IPL_DEPTH_32F,1);
IplImage* tmp_image = cvCreateImage( img_sz, IPL_DEPTH_32F,1 );
int corner_count =MAX_CORNERS;
CvPoint2D32f* cornersA = new CvPoint2D32f[MAX_CORNERS ];
cvGoodFeaturesToTrack(
imgA,
eig_image,
tmp_i.mage,
cornersA,
＆corner_count,
0.01,
5.0,
0,
3,
0,
0.04
）；
cvFindCornerSubPix(
imgA,
cornersA,
corner_count,
cvSize(win_size,win_size),
cvSize(-1,-1),
cvTermCriteria (CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,20,0.03)
）；
／／ Call Lucas Kanade algorithm
／／
char features_found[ MAX_CORNERS ];
float feature_errors[ MAX_CORNERS ];
CvSize pyr_sz = cvSize( imgA->width+8, imgB->height/3 );
IplImage* pyrA = cvCreateImage( pyr_sz, IPL_DEPTH_32F, 1 );
IplImage* pyrB = cvCreateImage( pyr_sz,IPL_DEPTH_32F,1);
CvPoint2D32f* cornersB = new CvPoint2D32f[ MAX_CORNERS ];
cvcalcoptica1FlowPyrLK(
imgA,
imgB,
pyrA,
pyrB,
cornersA,
CornersB,
corner_count,
cvSize( win_size,win_size ),
5,
features_found,
feature_errors,
cvTermCriteria( CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20,.3 ),
0
）；
／／ Now make some image of what we are looking at:
／／
for( int i=0; i<corner_count; i++){
if( features_found[i]==0|| feature_errors[i]>550 ){
printf("Error is &f/n", feature_errors [i]);
continue;
｛
printf("Got it/n");
CvPoint po = cvPoint(
cvRound( cornersA[i].x ),
cvRound( cornersA[i].Y )
）；
CvPoint p1 = cvPoint(
cvRound( cornersB[i].x),
CVRound( cornersB[i].Y)
）；
cvLine( imgC, p0, p1, CV_RGB(255,0,0),2 );
｛
cvNamedWindow("ImageA",0);
cvNamedWindow("ImageB",0);
cvNamedWindow("LKpyr_Opt icalFlow" , 0) ;
cvShowImage("ImageA" , imgA) ;
cvShowImage ( "ImageB" , imgB);
cvShowImage("LKpyr_OpticalFlow",imgC);
cvWaitKey(0);
return 0;
稠密跟踪方法
OpenCV 中还有另外两种光流方法，但现在已经很少用到了。这些函数比 Lucas-
Kanade 方法慢很多；另外，它们不支持图像金字塔匹配而不能用于跟踪大幅度的
运动。本节我们简单讨论一下这些函数。
【332~334】
Horn-Schunck方法
Horn 和 Schunck于1981年提出这种方法。此方法是首次使用亮度恒定假设和推导
出基本的亮度恒定方程的方法之一。Horn 和 Schunck 求解方程方法是假定一个速
度v,和v,的平滑约束。该约束是通过对光流速度分量的二阶导数进行规则化获
得：
0=(7+x7+x)27
0=(7+ay+a7)77-e@
这里α是不变的权重系数，称为规则化常数（regularization constant).α的值较大可
）
以获得更平滑（也就是更局部一致）的运动流向量。这是强制进行平滑的一个简单的
约束，其效果是惩罚光流变化剧烈的区域。与 Lucas-Kanade 算法一样，Horn-
Schunck 方法也要通过迭代来解微分方程。
void cvCalcopticalFloWHS(
const CvArr*
imgA,
const CvArr*
imgB,
int
usePrevious,
CvArr*
velx,
CvArr*
vely,
double
lambda,
CvTermCriteria
criteria
）；
【335~336】
其中 imgA和imgB必须是8位的单通道图像。velx和vely存放x和y方向的速
度，它们必须是32位的浮点数单通道图像。参数 usePrevious 指定算法使用从前
一帧计算的 velx和 vely 速度作为计算新的速度的初始值。参数 lambda 是与
Lagrange 乘子相关的权重。读者可能会问：“什么是Lagrange 乘子？”Lagrange
乘子出现在当我们（同时）最小化运动－亮度方程和平滑方程时；它表示当最小化时
赋予每一个方程误差的相对权重。
块匹配方法
读者可能会想：“光流到底有什么用？只是匹配上一帧和下一帧的像素而已。”有
一些人的确是这样用的。“块匹配”是对一类将图像分割成被称为“块”
［Huang95;Beauchemin95]的小区域的相似算法的统称。典型的块是正方形的，包
含了一定数目的像素。这些块可以重叠，实际中它们通常是重叠的。块匹配算法将
前一帧图像和当前图像划分成小块，然后计算这些块的运动。这一类算法在视频压
缩算法和计算机视觉的光流中扮演重要角色。
块匹配算法对像素的集合进行处理而非单个像素，所以返回的“速度图像”通常比
输入图像的分辨率低，但也有例外，这取决于块之间的重叠程度。下面的公式给出
了结果图像的尺寸：
①
这种情况下最好忽略这里的描述而设置1ambda为1.
W result=
Weo-Woec+
H resulr
Hoo-Haoe+Haee
ftoor
OpenCV 中的实现算法用螺旋搜索，找出（前一帧图像的）原始块的位置，然后再将
候选的新块同原始块进行比较。比较的结果是像素绝对差值的和（也就是L1距
离）找到很好的匹配后就终止搜索。函数原型如下：
【336】
id cvcalcopticalFlowBM(
const CvArr*
prev,
const CvArr*
curr,
CvSize
block_size,
CvSize
shift_size,
CvSize
max_range,
int
use_previous,
CvArr*
velx,
CvArr*
vely
）；
参数表示的意思简明易懂。参数prev和 curr 是前一帧图像和当前图像，它们都
是8位单通道的图像。block_size 是块的尺寸，shift_size 是块间移动的步长
（这个参数决定了块与块之间是否重叠以及如果重叠，重叠的程度又如何）。参数
max_range
是一个块周围邻域的尺寸，函数在下一帧中搜索此邻域以找到对应的
块。如果设置了参数use_previous,则 velx 和 vely 的值将作为块搜索的初始
值?.参数 velx和 vely是存储计算得的块的运动的32位的单通道图像。如前文
所说，运动的计算是在块级上进行的，所以输出结果的图像的坐标是对块而言的
（就是说，像素的集合体），而不是原始图像中的每个像素。
mean-shift和
mshift 跟踪
这一节中将要介绍两种跟踪方法，mean-shift 和 camshift(“camshaift”是
“continuously adaptive mean-shift”的缩写）。前者是可用于多种应用的通用的数据
①
如果use_previous == 0,块搜索区域与原始块的位置的距离为 max_range.如果
use_previous !=0,则块搜索区域的中心将被偏移Δx=vel(x,y)和
Δy=vel,(x,y).
分析方法（在第9章的分割问题中讨论过），计算机视觉正是这些应用之一。在介绍
mean-shift 的基本原理之后，会介绍 OpenCV 是如何用它在图像中进行跟踪运动
的。后一种方法camshift是建立在 mean-shift之上，它可以跟踪视频中尺寸可能产
生变化的目标。
mean-shift
mean-shift 算法①是一种在一组数据的密度分布中寻找局部极值的稳定的方法。若
分布是连续的，处理过程就比较容易，这种情况下本质上只需要对数据的密度直方
图应用爬山算法即可?.然而，对于离散的数据集，这个问题在某种程度上是比较
麻烦的。
【337】
这里用“稳定”这个形容词是指统计意义上的稳定，因为 mean-shift 忽略了数据中
的 outliers,即忽略远离数据峰值的点。mean-shift仅对数据局部窗口中的点进行处
理，处理完成后再移动窗口。
mean-shift算法的步骤如下。
（1)选择搜索窗口。
窗口的初始位置；
?窗口的类型（均匀、多项式、指数或者高斯类型）；
?窗口的形状（对称的或歪斜的，可能旋转的，圆形或矩形）；
?窗口的大小（超出窗口大小则被截去）。
（2)计算窗口（可能是带权重的）的重心。
（3)将窗口的中心设置在计算出的重心处。
（4)返回第（2)步，直到窗口的位置不再变化（通常会）?.
mean-shift 算法更加正式的描述是，mean-shift 算法与核密度估计的规则有关。
①
mean-shift是一个相当深奥的话题，我们在这里讨论的目的主要是使读者有个直观的印
象。参考Fukunaga[Fukunaga90]、Comaniciu 和 Meer[Comaniciu99]可以了解 mean-shift
的由来。
②
在这里使用“本质上”这个词是由于 mean-shift 会受尺度变化的影响。准确说就是：
mean-shift等价于先对连续分布用。mean-shift核进行卷积，然后再应用爬山算法。
③
迭代过程通常由最大迭代次数或者两次迭代中心变化的程度进行限制；虽然如此，迭
代过程最后都会收敛。
“核”是一个局部函数（例如高斯分布）。如果在足够的点处有足够合适的带权重和
尺度的核，数据的分布便可以完全根据这些核来表示。与核密度估计不同，mean-
shift 只估计数据分布的梯度（变化的方向）。若变化为0的地方则表示是这个分布的
峰值（虽然可能是局部的）。当然在附近或其他尺度上可能还有峰值。
图10-11是mean-shift算法中的计算公式。
xx
Start with akernel K(x-X)=ck
approximation of a probabllty
dsrulbuton P()=1/1
K(x-x).Focusor,the greadient VP(x)=/≥K(x-x).
Let:g(x)=-k'(x),the derivative of the kernel and we get:
VP(2)=K1/8
Meanshift
vector
图10-11:mean-shift等式及其含义
这些计算公式可用一个矩形的核进行简化，将 mean-shift 矢量等式简化为计算图
像像素分布的重心：
x.=M
y=Ma
这里，零阶矩的计算如下：
＝00W
ΣΣI(x,y)
一阶矩的计算为：
【338】
M0=ΣΣx(x,2)
和
Mo1=ΣΣyI(x,2)
①
矩形核并不随着到中心的距离下降，而是一个突然变成零的突然转换。这个与高斯核
的指数衰减不同，与Epanechnikov核的随着到中心的距离的开方衰减也不同。
mean-shift 矢量告诉我们如何将 mean-shift 窗口的中心重新移动到由计算得出的此
窗口的重心的位置。很显然，窗口的移动造成了窗口内容的改变，于是我们又重复
刚才重新定位窗口中心的步骤。窗口中心重定位的过程通常会收敛到mean-shift矢
量为0(也就是，窗口不能再移动）。收敛的位置在窗口中像素分布的局部最大值（峰
值）处。由于“峰值”本身是一个对尺度变化敏感的量，所以窗口大小不同，峰值
的位置也不一样。
图 10-12 是一个二维数据分布及初始化窗口（此例中，窗口为矩形）的例子。箭头表
示迭代过程最终收敛到分布的局部峰值。正如前文所说的，我们可以看到由于处于
mean-shift 窗口外部点不影响 mean-shift的收敛性，所以寻找峰值的这一过程在统
计意义上是稳定的。
图10-12:mean-shift算法的过程。将初始窗口放置在数据点的二维数组上，
连续重定位（窗口）中心到数据分布的局部峰值处，直到收敛
1998年，人们认识到这种寻优的算法可以用于跟踪视频中的运动物体［Bradski98a;
Bradski98b],此后这个算法被大大地扩展［Comaniciu03].OpenCV 中运用 mean-
shift 算法的函数是在图像分析的层面上实现的。这也就是说，OpenCV 中实现的
mean-shift 把一幅代表将要被分析的密度分布的图像作为输入，而不是一组任意的
点（可能是任意的维数）作为输入。
【339~340】
int cvMeanShift(
const CvArr*
prob_image,
存在mhi中的最大持续时间。换句话说，mhi中任何比timestamp 减去 duration
的值早（少）的像素将被置为0.
一旦运动模板记录了不同时间的物体轮廓，就可以用计算 mhi 图像的梯度来获取
全局运动信息。计算出的这些梯度（例如，用第6章讨论过的 Scharr 或 Sobel 梯度
函数来计算），一些梯度值会很大并且是无效的。mhi 图像中被置0的旧的或没有
运动的部分的梯度是无效的，因为它们会在轮廓的外部边缘区域人为地产生很大的
梯度值，如图10-15(A)所示。由于使用cvUpdateMotionHistory()将新的轮廓引
进 mhi的时间步长是已知的，所以梯度值应该为多大也就知道了（即为 dx和 dy
step derivatives).那么我们可以用梯度幅值来消除特别大的梯度，如图10-15(B)所
示。最终，我们就得到全局运动的测量，如图 10-15(C).图 10-15(A)和（B)是由
cvCalcMotionGradient()函数得到的：
【343】
void cvCalcMotionGradient(
const CvArr*
mhi,
CvArr*
mask,
CvArr*
orientation,
double
deltal,
double
delta2,
int
aperture_size=3
）；
（B)
C
图10-15:mhi图像的运动梯度：（A)梯度值和方向；（B)大的梯度被去除；（C)
得到全局运动方向
函数cvCalcMotionGradient 中，所有的图像矩阵都是单通道的。函数的输入
mhi是一幅浮点数值的运动历史图像，delta1和delta2则分别是允许的最小和最
大梯度值。期望的梯度值是连续调用 cvUpdateMotionHistory()的每个轮廓间当
前时间标记的平均数；将delta1 设置为小于此平均值的一半，delta2 设置为大
于此平均值的一半比较合适。变量 aperture_size 设置梯度算子的宽和高，其值
可以为－1(3x3 CV_SCHARR梯度滤波器）、3(默认的3x3 Sobel滤波器）、5(5x5 Sobel
滤波器）或者7(7x7滤波器）。函数的输出mask是一幅8位的单通道图像，其中的
378 第10章
CvRect
window,
CvTermCriteria
criteria,
CvConnectedComp*
comp
）；
CVMeanShift()函数中，参数 prob_image是单通道图像，代表可能位置的密度，
其类型可以是 byte 或 float.参数 window的位置设置在为指定初始位置，大小为
核窗口的大小。终止条件criteria 主要由 mean-shift移动的最大迭代次数和可视
为窗口位置收敛的最小移动距离组成，其定义在其他章节中有描述①.连接部件
comp 在comp->rect中包含收敛后搜索窗口的位置，在comp->area中包含了窗口
内部所有像素点的和。
函数cvMeanShift()使用的是矩形窗口的mean-shift 算法，但这同样可以用于跟
踪。在这种情况下，需要首先选择代表物体的特征的分布（例如，颜色＋纹理），然
后在物体的特征分布上开始mean-shift窗口搜索，最后计算下一帧视频中所选择的
特征的分布。从当前的窗口位置开始，mean-shift 算法寻找特征分布的新的峰值或
者 mode,它们（被假定）设置在最初产生颜色和纹理的物体的中心。这样，mean-
shift窗口就可以逐帧跟踪物体的运动。
CamShift
另一个相关的算法是CamShift跟踪器。与mean-shift 不同的是，CamShift 搜索窗
口会自我调整尺寸。如果有一个易于分割的分布（例如保持紧密的人脸特征），此算
法可以根据人在走近或远离摄像机时脸的尺寸而自动调整窗口的尺寸。CamShift
算法的形式如下：
?int cvCamShift(
const CvArr*
prob_image,
CVRect
window,
CvTermCriteria
criteria,
CvConnectedComp* comp,
CVBox2D*
box = NULL
）；
前四个参数的含义与cVMeanShift()算法中的含义相同。如果参数box 不为空，
①
强调一下mean-shift 总是会收敛，但是如果一个分布在局部峰值附近非常“平坦”，
收敛将会变得非常慢。
Current
TOIE
SIllouette
Current
Current
Sllouette
Sillouette
Current
则它包含了新的尺寸的box,其中也包含根据二阶矩计算出的物体的方向。在跟踪
的应用中，会把由前一帧计算出的新尺寸的box作为下一帧的window.
注意：很多人认为mean-shift 和 camshift 利用颜色特征进行跟踪，但实际上不完全是这
样。两种算法跟踪在prob_image 中所代表的任意的分布，所以它们是轻量级的、稳
定和有效的跟踪器。
运动模板
运动模板（motion template)是由 MIT 媒体实验室的 Bobick 和 Davis 提出的
［Bobick96;Davis97],并且其中的一人又参与了它进一步的开发［Davis99;
Bradski00].而这一步工作为运动模板在OpenCV 中的实现提供了基础。运动模板
是一种有效的跟踪普通运动的方法，尤其可应用在姿态识别中。运用运动模板需要
知道物体的轮廓（或者轮廓的一部分，此处轮廓指 silhouette,即实心轮廓），而轮廓
的获取有不同的方法。
【341~342】
（1)最简单的获取物体轮廓的方法是用一个静止的摄像机，再使用帧间差（见第9
章论述）得到物体的运动边缘，这种信息就足够让运动模板发挥效用。
（2).利用颜色信息。例如，如果已知背景的颜色，比如亮绿色，那么可以很简单地
将不是亮绿色的任何物体看作前景。
（3)另一种方法（在第9章中讨论过）是学习一个背景模型，从此背景中可以将新的
前景物体／人以轮廓的形式分割出来。
（4)使用主动轮廓技术。例如，创建一个近红外光的墙，将能感应近红外线的摄像
机对着墙，那么任何介于墙与摄像机之间的物体都会以轮廓的形式出现。
（5)使用热感图像。任何温度高的物体（如人脸）都可以被认为是前景。
（6)最后，可以用第9章中谈到的分割技术（例如金字塔分割或 mean-shift分割）来
生成轮廓。
从现在开始，假设我们有一个如图10-13(A)中白色矩形所代表的被较好的分割出的
物体轮廓。白色表示所有像素都被设置为最新的系统时间的浮点数值。随着矩形的
运动，新的轮廓将被捕获并且被（新的）当前轮廓覆盖；新的轮廓为图 10-13(B)和
图 10-13(C)中的白色矩形所示。较早的运动在图10-13中表示为亮度渐暗的矩形。
这些连续变暗的轮廓记录了早前运动的历史，所以被称为“运动历史图像（motion
history image)”。
图10-13:运动模板图：（A)当前时间分割出的物体（白色）；（B)在下一个时间
点，物体运动并以（新的）当前时间标记，将前面的分割边界留在后面；（C)在下
一个时间点，物体继续运动并将此前的分割标记为渐暗的矩形，这些包含运动
的序列就产生运动历史图像
【342】
若轮廓的时间超过设定的比当前时间早的持续时间，则其被置0,如图10-14所
示。OpenCV中完成运动模板构建的函数是cvUpdateMotionHistory():
void cvupdateMotionHistory(
const CvArr*
silhouette,
CvArr*
mhi,
double
timestamp,
double
duration
）；
（A)
（B)
运动
运动
图10-14:两个物体的运动模板轮廓（A);超过设定持续时间的轮廓被置0(B)
在cvupdateMotionHistory()中，所有图像都是单通道图像。图像
silhouette
是一幅单字节图像，其中非0像素代表前景物体的最新的分割轮廓。图像 mhi是
一幅浮点值的图像，代表运动模板（也就是运动历史图像）。timestamp 是当前系统
时间（典型地以毫秒为单位），duration 如刚才谈到的表示运动历史像素的允许保
第10章
跟踪与运动
跟踪基础
当我们在处理一段视频而非某张静止的图像时，通常视频中有一个或几个特定物体
是我们想在整个视野范围内关注的。在前一章中，我们知道了如何在逐帧（frame-
by-frame)分离出特定的形状，例如人或者车；现在，我们尝试理解这个物体的运
动。理解物体的运动则主要包含两个部分：识别（identification)和建模。
识别指在视频流后续的帧中找出之前某帧中的感兴趣物体。前面章节中讲到的矩和
颜色直方图可以帮助识别我们关注的物体。一个相关的问题就是跟踪不明物体。当
我们需要根据运动来确定什么是感兴趣的，或者当物体的运动使得其成为感兴趣物
体时，跟踪不明物体就很重要了。经典的跟踪不明物体的方法是跟踪视觉上重要的
关键点，并不是整个物体。OpenCV 提供了两种方法实现跟踪关键点：Lucas-
Kanade?[Lucas81]和 Horn-Schunk [Horn81]方法。这两种方法分别代表了通常提到
的稀疏（sparse)和稠密（dense)光流。
上述方法实际上只能给出物体实际位置的初步计算。理解物体运动的第二个部分，
也就是建模，则可以帮助我们解决这个问题。为估计由这种粗略测量得出的物体的
轨迹，人们设计了许多有力的数学方法。这些方法适用于二维和三维物体或物体位
置的模型。
①很奇怪，OpenCV 中实现的金字塔框架下的Lucas-Kanade 光流方法是根据一篇未正式
发表的论文实现的，这篇论文是Bouguet的［Bouguet04]:
寻找角点
有多种局部特征可以用来进行跟踪。花一点时间去思考究竟角点（corner)用什么样
的特征来描述是值得的。很明显，如果从一面很大的空墙上选择一个点，那么将很
难从视频的下一帧中再找出这一点。如果墙上的所有点都是一样的或者是相似的，
我们就不会有太好的运气能在随后的视频帧中跟踪到这个点了。相反，如果选择一
个独一无二的点，那么再找到这点的几率就非常大。实际上，被选择的点或特征应
该是独一无二的，或者至少接近独一无二，并且在与另一张图像的其他点可以进行
参数化的比较。如图10-1所示。
【316~317】
图10-1:圆圈圈出的特征点是易于跟踪的点，而矩形圈出的特征－
甚至是明
显的边缘－却不然
回到那面大而空的墙，我们想要找本身具有明显变化的点－例如导数值比较明显
的地方。虽然仅此是不够的，但这是一个开始。一个导数值比较明显的点可能在某
种类型的缘上，但却和此边缘上的其他点看起来一样（参考图 10-8和“Lucas-
Kanade方法”一节中讨论的孔径问题）。
如果一个点在两个正交的方向上都有明显的导数，则我们认为此点更倾向是独一无
二的，所以，许多可跟踪的特征点都称为角点。从直观上讲，角点（而非边缘）是一
类含有足够信息且能从当前帧和下一帧中都能提取出来的点。
最普遍使用的角点定义是由Harris[Harris88]提出的。定义的基础是图像灰度强度
的二阶导数（x,?2y,dxdy)矩阵。考虑到图像所有的像素点，我们可以想像图像的二
阶导数即形成一幅新的“二阶导数图像”，或融合在一起形成一幅新的 Hessian 图
像。这个术语源自于一个点的如下定义的二维Hessian矩阵：
跟踪与运动
H(p)=
【317】
对于 Harris 角点，我们使用每点周围小窗口的二阶导数图像的自相关矩阵。这个
自相关矩阵的定义如下：
2(x+iy+j)
（x+i,y+j),(x+y+17
M(x,y)=
3
w.1(x+i,y+j)I,(x+i,y+j)
（x+3y+)
－KSI,jsK
这里w,是可以归一化的权重比例，但是通常被用作产生圆形窗口或高斯权重。
Harris 定义的角点位于图像二阶导数的自相关矩阵有两个最大特征值的地方，这在
本质上表示以此点为中心周围存在至少两个不同方向的纹理（或者边缘），正如实际
的角点是由至少两个边缘相交于一点而产生。之所以采用二阶导数是由于它对均匀
梯度不产生响应?.角点的这个定义还有另一个优点。被跟踪的物体在移动过程中
也可能会旋转，找到同时对移动和旋转不变的量是很重要的。只考虑自相关矩阵的
特征值可达到这个目的。这两个最大特征值不仅可以判定一个点是否一个好的可跟
踪的特征点，同时也提供了对这个点进行识别的一个标识。
Harris 最原始的定义是将矩阵H(p)的行列式值与H(p)的迹（带权重系数）相减，再
将差值同预先给定的阈值进行比较。后来 Shi 和Tomasi[Shi94]发现，若两个特征
值中较小的一个大于最小阈值，则会得到强角点。Shi 和 Tomasi 的方法比较充
分，并且在很多情况可以得到比使用Harris方法更好的结果。
函数cvGoodFeaturesToTrack()采用 Shi和 Tomasi提出的方法，先计算二阶导数
（利用Sobel算子），再计算特征值，它返回满足易于跟踪的定义的一系列点。
void cvGoodFeaturesToTrack(
const CvArr*
image,
CvArr*
eigImage,
CvArr*
tempImage,
CvPoint2D32f*
corners,
int*
corner_count,
double
quality_level,
double
min_distance,
①
梯度是一阶导数而来。若一阶导数是均匀的（常数），则二阶导数为0.
st
CvArr*
mask
＝NULL,
block_size
＝3,
use_harris
＝0,
ouble
k
＝0.4
【318】
数中，
输入图像 image须为8位或32位（也就是，IPL_DEPTH_8U 或者
32F)单通道图像。第二和第三个参数是大小与输人图像相同的32位单
参数
tempImage 和 eigImage 在计算过程中被当作临时变量使用，计
量eigImage中的内容是有效的。特别地，每个元素包含了输入图像中
最小特征值。
corners 是函数的输出，为检测到的32位（CvPoint2D32f)
组，在调用函数cvGoodFeaturesTorrack()之前需要给这个数组分配内
很明显，只能为此数组分配有限的内存。corner_count表示可以返回的
数目。函数调用结束后，corner_count 输出实际检测到的角点数目。
level
表示一点被认为是角点的可接受的最小特征值。实际用于过滤角
最小特征值是
quality_level 与图像中最大特征值的乘积。所以
1eve1的值不应超过1(常用的值为0.10或者0.01).检测完所有的角点
要进一步剔除掉一些距离较近的角点。min_distance 保证返回的角点之间
5不小于min_distance个像素。
是可选参数，是一幅像素值为布尔类型的图像，用于指定输人图像中参与角
的像素点。若mask的值设为NOLL,则选择整个图像。block_size是计算
自相关矩阵时指定点的领域，采用小窗口计算的结果比单点（也就是
size为1)计算的结果要好。若use_harris的值为非0,则函数使用
的角点定义；若为0,则使用shi-Tomasi的定义。当 use_harris为k且
则k为用于设置Hessian 自相关矩阵即对 Hessian行列式的相对权重的权重
cvGoodreaturesTorrack()的输出结果为角点的位置数组，我们往往会希望
第一幅相似的图像中寻找这些角点。就本章而言，我们关注的是在后续的视频帧
我到这些特征点，但是这个函数也有其他的应用，例如，对从微小不同的视角拍
多幅图像进行匹配。这个问题将在后面讨论立体视觉的小节中遇到。
花
亚像素级角点
如果我们进行图像处理的目的不是提取用于识别的特征点而是进行几何测量，则通
常需要更高的精度，而函数cVGoodFeaturesToTrack()只能提供简单的像素的坐
标值，就是说，我们有时会需要实数坐标值而不是整数坐标值，例如（8.25,
117.16).
您也许需要确定图像中一个尖锐的峰值点的位置，但是峰值的位置一般都不会恰好
位于一个像素的正中心，使得确定其位置比较困难。要解决这个问题，需要先用一
条曲线（例如一条抛物线）拟合图像的值，再用数学的方法计算位于像素之间的峰值
点的位置。亚像素检测方法就是一些有关峰值点位置的计算技巧（了解更多此方面
知识以及最新的方法，请参考Lucchese[Lucchese02]和Chen[Chen05]).图像测
量常用的领域为三维重建、摄像机标定、图像拼接以及在卫星图像中查找特定信
号，如一栋建筑的精确位置。
【319~320】
亚像素级角点的位置在摄像机标定、跟踪并重建摄像机的轨迹或者重建被跟踪目标
的三维结构时是一个基本的测量值。我们已经知道如何对整数网格像素来检测角点
并得到角点的位置，下面所要讨论的是如何将求得的角点位置精确到亚像素级精
度。一个向量和与其正交的向量的点积为0,角点则满足这种情况，如图10-2
所示。
邻域
边上的梯度
VI(P)
0=(d)I
（a):p点在一个区域内部
（b):p点在区域边缘
在这两种情况下，点积都为零
＜VI(p),q-p>=0
图10-2:计算亚像素级精度的角点：（a)点p附近的图像是均匀的，其梯度为
0;(b)边缘的梯度与沿边缘方向的q-p向量正交。在图中两种情况下，p点的
梯度与q-p向量的点积均为0(参考正文）
图10-2假设起始角点q在实际亚像素级角点的附近。检测所有的q-p向量。若点p
位于一个均匀的区域，则点p处的梯度为0.若q-p 向量的方向与边缘的方向一
致，则此边缘上p点处的梯度与q-p向量正交，在这两种情况下，p点处的梯度与
q-p向量的点积为0.我们可以在p点周围找到很多组梯度以及相关的向量q-p,令
其点集为0,然后可以通过求解方程组，方程组的解即为角点q的亚像素级精度的
位置，也就是精确的角点位置。
【320】
函数cvFindCornerSubPix()用于发现亚像素精度的角点位置。
void cvFindCornerSubPix(
const CvArr*
image,
CvPoint2D32f*
corners,
int
count,
CvSize
win,
CvSize
zero_zone,
CvTermCriteria
criteria
）；
输入图像 image是8位单通道的灰度图像。corners 为整数值的像素位置，例如
是由cVGoodFeatures()得到的位置坐标。corners
设定了角点的初始位置。
count为需要计算的角点数目。
实际计算亚像素级的角点位置时，解的是一个点积的表达式为0的方程组（参考图
10-2),其中每一个方程都是由q邻域的一个点产生。win 指定了等式产生的窗口
的尺寸。搜索窗口的中心是整数坐标值的角点并从中心点在每个方向上扩展 win
指定的像素（也就是说，如果win.width=4,则实际的窗口大小为4+1+4=9个像
素宽）。这些等式构成一个可用自相关矩阵的逆（非前文讨论 Harris 角点时涉及的
自相关矩阵）来求解的线形方程组。实际上，由于非常接近p的像素产生了很小的
特征值，所以这个自相关矩阵并不总是可逆的。为了解决这个问题，一般可以简单
地剔除离p点非常近的像素。输入参数zero_zone定义了一个禁区（与win相似，
但通常比 win小），这个区域在方程组以及自相关矩阵中不被考虑。如果不需要这
样一个禁区，则zero_zone应设置为cvSize(-1,-1).
当找到一个q的新位置时，算法会以这个新的角点作为初始点进行迭代直到满足用
户定义的迭代终止条件。迭代过程的终止条件可以是最大迭代次数
CV_TERMCRIT_ITER 类型，或者是设定的精度 CV_TERMCRIT_ITER 类型（或者是两
者的组合）。终止条件的设置在极大程度上影响最终得到的亚像素值的精度。例
如，如果指定0.10,则求得的亚像素级精度为像素的十分之一。
不变特征
继 Harris 提出角点以及后来 Shi和 Tomasi提出角点后，许多其他类型的角点和相
关局部特征点也被提出来。SIFT(“scale-invariant feature transform”）是其中一种
广泛应用的类型［Lowe04].SIFT 特征正如其名称一样是缩放不变的。SIFT 在一点
处检测主要梯度方向，根据这个方向记录局部梯度直方图结果，所以 SIFT 也是旋
转不变的。正是由于这些特点，SIFT 特征在小的仿射变换中有相对不错的表现。
虽然OpenCV库中还没有实现SIFT算法（参考第14章），但可以用OpenCV的基本
函数实现它。我们将不在这个问题上花更多的时间，但需要记住的是，用前面已经
介绍过的 OpenCV 函数可以实现计算机视觉领域的文献中提出的绝大部分特征。
【321~322】
光流
如前所述，我们经常不知道任何关于视频内容的先验知识，但是需要估计两帧之间
（或一个帧序列）的运动。运动本身就说明了存在运动的感兴趣的目标。在图10-3
中解释了光流。
a
光流
2
23
I(t+1)
（a)/(t),点集
：｛pi}
（b)速度向量｛3}
图像跟踪
（c)图像序列（单个摄像机）
（d)跟踪结果
图10-3:光流。目标特征点（a)被连续跟踪，它们的运动被转换成速度矢量
（b);下面的两幅图分别表示ー幅走廊的图像（c)和摄像机沿着走廊运动时的光
流矢量（d)(原图由 Jean-Yves Bouguet提供）
可以将图像中的每个像素与速度关联，或者等价地，与表示像素在连续两帧之间的
位移关联。这样得到的是稠密光流（dense optical flow),即将图像中的每个像素都
数
与速度关联。Horm-Schunck 方法［Horn81]计算的就是稠密光流的速度场。OpenCV
中实现了一种比较简单直接的方法，即对前后连续两帧的一个像素的邻域进行匹
配。这种方法被称为块匹配（block matching).这两种方法都会在“稠密跟踪方法”
一节中讨论。
【322~323】
实际上计算稠密光流并不容易。我们先来看看一张白纸的运动。上一帧中的白色的
像素在下一帧中仍然为色，而只有边缘的像素并且是与运动方向垂直的像素才会
产生变化。稠密光流的方法需要使用某种插值方法在比较容易跟踪的像素之间进行
插值以解决那些运动不明确的像素，从中可以清楚地看到稠密光流相当大的计
算量。
于是我们找到了另一种方法，稀疏光流（sparse optical flow).稀疏光流的计算需要
在被跟踪之前指定一组点。如果这些点具有某种明显的特性，例如前面讲到的“角
点”，那么跟踪就会相对稳定和可靠。OpenCV 可以帮助我们找到最适合跟踪的特
征点。在很多的实际应用中，稀疏跟踪的计算开销比稠密跟踪小得多，以至于后者
只有在理论中研究。
下面的几节介绍了几种不同的跟踪的方法。我们从最流行的稀疏跟踪方法 Lucas-
Kanade(LK)光流开始介绍；这种方法与图像金字塔一起，可以跟踪更快的运动。
接下来介绍两种稠密跟踪方法，Horn-Schunck方法和块匹配方法。
Lucas-Kanade 方法
Lucas-Kanade 算法（LK)[Lucas81]是用于求稠密光流的，最初于1981年提出。由于
算法易于应用在输入图像中的一组点上，其后来成为求稀疏光流的一种重要方法。
LK 算法只需要每个感兴趣点周围小窗口的局部信息，所以它可以应用于稀疏内
容，这与 Horn 和 Schunck 的算法全局性是不同的（后面会详细讲述）。但是，使用
小窗口的LK算法存在不足之处，较大的运动会将点移出这个小窗口，从而造成算
法无法再找到这些点。金字塔的LK算法可以解决这个问题，即从图像金字塔的最
高层（细节最少）开始向金字塔的低层（丰富的细节）进行跟踪。跟踪图像金字塔允许
小窗口捕获较大的运动。
这个算法是一个重要而有效的方法，我们将会具体介绍其中的数学问题。不想阅读
①
Black 和 Anadan 提出的计算稠密光流的方法［Black93;Black96]被广泛用在电影制作
中，因为为了追求视觉的质量，电影工作室愿意牺牲时间开销来获得详细的光流信
息。OpenCV未来的版本中会将这些方法包含进去（参考第14章）。
这些细节的读者可以直接阅读函数描述和代码部分，但是建议读者至少浏览一遍其
中描述Lucas-Kanade 光流假设条件的文字和插图部分，这样在失去头绪的时候知
道应该如何去解决。
Lucas-Kanade 算法原理
LK算法基于以下三个假设。
（1)亮度恒定。图像场景中目标的像素在帧间运动时外观上保持不变。对于灰度图
像（LK算法也可用于彩色图像），需要假设像素被逐帧跟踪时其亮度不发生变化。
（2)时间连续或者运动是“小运动”。图像的运动随时间的变化比较缓慢。实际应
用中指的是时间变化相对图像中运动的比例要足够小，这样目标在帧间的运动
就比较小。
（3)空间一致。一个场景中同一表面上邻近的点具有相似的运动，在图像平面上的
投影也在邻近区域。
图10-4 说明了这三点假设如何构建一个有效的跟踪算法。第一点假设：亮度恒
定，指对被跟踪部分（patch)像素不随时间变化的要求：
（p+1*(p+1)x)I=(r(x)x)I=(r*x)f
假设
亮度恒定
（a)
1(x+u,y+v,t+1)=1(x,y,t)
时间连续
空间一致
表面
（q)
（c)
图像平面
图10-4:Lucas-Kanade 光流假设：（a)场景中物体被跟踪的部分的亮度不变；
（b)运动相对于帧率是缓慢的；（c)[相邻的点保持相邻（图由 Michael Black
［Black82]提供）］
【324】
很简单，公式的意思是被跟踪像素的灰度不随时间变化而变化：
af(x)=0
10
第二点假设：时间连续，指相邻帧之间的运动较小。换句话，可以将运动的变化看
成是亮度对时间的导数（即断言一个序列中相邻帧间的变化是小微分的）。我们先来
看一下一维空间的例子，以理解第二点假设所蕴含的意思。
由亮度恒定的等式开始，考虑隐含的x为t的函数这点，将亮度的定义f(x,t)用
I(x(t),t)替换，再应用偏微分的链式规则（chain rule),即：
o-
其中，I,是图像的偏导数，I,是图像随时间的导数，v是要求的速度。简单一维空
间中的光流速度的等式为：
＝－1/1
下面我们对一维跟踪问题进行一些直观的解释。图 10-5 表示的是左边值大右边值
小、沿x轴向右运动的一个“边缘”。我们要求出图的上半部分标出的边缘运动的
速度V.在图的下半部分，对这个速度的测量正是“rise over run”，即 rise 是随
时间变化的而run是斜率（空间的导数）。负号纠正了x的斜率。
一维光流
1(x,t)
［（x,t+1)
I(x,t)
1(x,t+1)
时域导数
空域导数
假定：
?亮度不变
／e=
＊小运动
图10-5:一维空间中的Lucas-Kanade光流。通过计算亮度对时域导数与亮度
对空域导数的比值，可以估计出运动边缘（图的上半部分）的速度
图10-5 揭示了光流公式的另一个问题：我们的假设可能不是十分正确，即图像的
亮度实际上不是恒定的，时间步长（由摄像机设定）也不是如我们期望的相对于运动
足够短。所以我们求解的速度并不准确。但是，如果求解到的速度与实际的速度
“足够接近”，就可以用迭代的方法来解决这个问题，如图 10-6所示，将第一次
估计的速度（不准确的）作为初始值进行下一次迭代并重复这个过程。根据沿x运动
的像素不变的亮度恒定假设，在迭代过程中，可以保持使用由第一帧计算得到的x
的空间导数。这种重复使用已经计算出的空间导数极大地节省了计算开销。时间导
数仍需要在每次迭代和每一帧中重新计算，但如果初始值与实际值足够接近，迭代
过程在5次迭代之内就会收敛到接近的准确值，这就是所谓的牛顿法。如果初始估
计值与实际值相差比较大，则牛顿法实际上是发散的。
【325~326】
代可以优化速度向量
I(x,t)
I(x,t+1)
第2次迭代得到的时域导数
x
J
空域导数也采用
Fpmerion 4/4
约5次迭代后收敛
图10-6:通过迭代来优化光流的求解（牛顿法）。使用同样的两个图像和同样的
空域导数，重新求解时域导数，通常可通过几次迭代收敛到稳定解
前面是LK算法在一维空间的用法，下面我们将其扩展应用到二维空间的图像上。
乍看起来好像比较简单，只需要加入y坐标即可。改变一下符号，速度的y分量为
v,x分量为u,则得到：
Iu+I,v+I,=0
不幸的是，这个等式对任意一个像素都含有两个未知量。这就是说对于单个像素，
等式的约束条件过少，不能得到此点二维运动的定解。我们只能求得与光流方程描
述的方向垂直的运动分量。图10-7给出了数学和几何的细节。
垂直光流由孔径问题产生，即用小孔或小窗口去测量运动。这种情况下，我们通常
只能观测到边缘而观测不到角点，而只依靠边缘是不足以判断整个物体是如何运动
（也就是朝哪个方向运动）的，如图10-8所示。
6九
从一维到二维光流
公式
但是，只能求得一条直线，而不是一个点
Ixu+Iyv+I4=0
Ixu+Iyv=-14
I-=n IΔ
1,u+1,+1=0
＝1A=
＂Normal flow”
图10-7:单个像素的二维光流：单个像素的光流是无法求解的，最多只能求出
光流方向，因为光流方向与光流方程描述的线垂直（图由Michael Black提供）
（a)
孔经问题
（b)
图10-8:孔径问题：从aperture window(a)我们可以观测到边缘向右运动，但
是无法观察到边缘也在向下运动（b)
那么，如何解决单个像素不能求解整个运动的问题呢？这时需要利用光流的最后一
点假设。若一个局部区域的像素运动是一致的，则可以建立邻域像素的系统方程来
求解中心像素的运动。例如，如果用当前像素 5x5①邻域的像素的亮度值（彩色像素
的光流只需增加两倍）来计算此像素的运动，则可以建立如下的25个方程。
【326~327】
①
窗口可以是 3x3、7x7 或者任何被指定的值。若窗口太大则会由于违背运动一致的假
设而不能进行较好的跟踪。若窗口太小，则又会产生孔径问题。
（
1,(p1)
I,(P2) I,(P2)
I,(P2s)
1,(P25)
25x2
现在我们得到一个约束条件过多的系统方程，若在5x5的窗口中包含两条或院上
边缘则可以解此系统方程。为了解这个系统方程，需要建立一个该方程的最小事
方，通过下面方程来求解最小化的｜Ad-bP:
／＝（
由这个关系式可以得到u和v运动分量。这个关系更详细的表述如下：
S11, ΣI1,u
ΣIIy
Σ1,
当（AA)可逆时，方程的解如下：
＝（AA)-1ATb
【328】
当（AA)满秩（秩为2)也即（AA)有两个较大特征向量时，（AA)可逆。图像中纹理至
少有两个方向的区域，这个条件可以满足。这种情况下，跟踪窗口的中心在图像的
角点区域时，（AA)的特性最好。由这一点我们可以联系到前面讨论的 Harris 角点
检测器。事实上，正因为（AA)在角点出有两个大的特征向量，所以这些角点是
“可用于跟踪的良好特征点”（参考前面对 cvGoodFeaturesToTrack()的介绍），
后文即将介绍的cvCalcopticalFlowLK()实现了这个计算。
理解了小而连贯运动的假设的读者现在会感到不解：对于大多数30Hz的接像机。
大而不连贯的运动是普遍存在的情况。而 Lucas-Kanade 光流正因为这个原因在实
际中的跟踪效果并不是很好：我们需要一个大的窗口来捕获大的运动，面大窗口往
往会违背运动连贯的假设！图像金字塔可以解决这个问题，即最初有铰大的空间尺
度上进行跟踪，再通过对图像金字塔向下直到图像像素的处理来修正初始运动速度
的假定。
所以，建议的跟踪方法是：在图像金字塔的最高层计算光流，用得到的运动基计结
果作为下一层金字塔的起始点，重复这个过程直到到达金字塔的最底层。这样就将
不满足运动假设的可能性降到最小从而实现对更快和更长的运动的跟踪。这个更精
致的算法叫金字塔Lucas-Kanade 光流，其原理在图10-9中描述。实现金字塔
Lucas-Kanade 光流的函数是 cvCalcOpticalFlowPyrLK(),后面将介绍这个函数。
从粗到精的光流估计
run iterative L-K
warp and upsample
run iterative L-K<
image1
高斯金字塔图像
Ic-1
高斯金字塔图像
图 10-9:金字塔Lucas-Kanade光流。为了减轻由小而连贯的运动假设引起的
问题，首先在金字塔顶层计算光流；在上一次估计到的运动将作为下一层的起
始点，进行进一步估计
Lucas-Kanade 代码
下面的函数实现了非金字塔的Lucas-Kanade稠密光流算法：
void cvCalcopticalFlowLK(.
const CvArr*
imgA,
const CvArr*
imgB,
CvSize
winSize,
CvArr*
velx,
CvArr*
vely
）；
OpenCV 这个函数的输出只记录了可以计算最小误差的像素。对于最小误差不能被
可靠地计算出的像素，相关的速度被设置为0.多数情况下我们不会使用这个函
数，而用下面介绍的基于图像金字塔的方法。
金字塔Lucas-Kanande 代码
现在来介绍 OpenCV 中在图像金字塔中计算 Lucas-Kanade 光流的算法
cvCalcopticalFlowPyrLK().我们将会看到，这个计算光流的函数使用了“易于
跟踪的特征点”并返回每个点被跟踪的情况。
【329~330】
void cvCalcopticalFlowPyrLK(
const CvArr*
imgA,
const CvArr*
imgB,
CvArr*
pyrA,
CvArr*
pyrB,
CvPoint2D32f*
featuresA,
CvPoint2D32f*
featuresB,
int
count,
CvSize
winSize,
int
level,
char*
status,
float*
track_error,
CvTermCriteria criteria,
int
flags
）；
此函数有很多输入参数，我们现花一点时间来弄清楚它们是什么。一旦掌握了这个
函数，我们就可以到下一个问题，即跟踪哪些点和怎样计算。
cvCalcopticalFlowPyrLk()的前两个参数代表初始图像和最终图像，两幅图像都
是8位的单通道图像。第三、四个参数是申请存放两幅输入图像（pyrA和 pyrB)金
字塔的缓存。这两个缓存的大小至少为（img.width+8)*img.height/3个字节?.
（如果这两个指针设置为空，则当函数被调用时会自动分配合适的内存，使用并释
放之，但这样会降低函数的执行效率。）featuresA 数组存放的是用于寻找运动的
点，featuresB 与 featuresB 相似，存放 featuresA 中点的新位置；count是
featuresA 中点的数目。win_size定义了计算局部连续运动的窗口尺寸。参数
level 用于设置构建的图像金字塔的栈的层数。若 level 设置为0则不使用金字
塔。数组 status的长度是 count;函数调用结束时，status 中的每个元素被置
1(对应点在第二幅图像中被发现）或0(对应点在第二幅图像中未被发现）。参数
track_error 为可选参数，它是表示被跟踪点的原始图像小区域与此点在第二幅
图像的小区域间的差的数组。track_error 可以删除那些局部外观小区域随点的
①分配这个大小的缓存是因为演算空间不仅存放图像本身，还要存放整个金字塔。
运动变化剧烈的点。
我们还需要一个迭代终止的条件，即 criteria.CvTermCriteria 是被许多含有
迭代过程的OpenCV算法使用的一个结构：
cvTermCriteria(
int type, // CV_TERMCRIT_ITER, CV_TERMCRIT_EPS, or both
int max_iter,
double epsilon
）；
我们使用cvTermCriteria()函数来生成需要的这个结构。函数的第一个参数为
CV_TERMCRIT_ITER 或者 CV_TERMCRIT_EPS.这个参数分别告诉算法是根据迭代
次数来终止迭代过程，还是根据收敛误差是否达到一个较小值来终止迭代。后面的
两个参数则设置终止算法的一个或两个准则的值。有两个选项的原因使我们可以设
税市桑CERCRIITERICV_TERMCRITEPS,刻兼市
时迭代终止（大多数的实际代码就是这样设置的）。
最后一个参数 flags允许对函数内部 bookkeeping进行一些细微的控制；它可以设
置为下面的任何一个或全部（使用位OR运算）。
CV_LKFLOW_PYR_A_READY 在调用和存储到
pyrA 之前，先计算第一帧
的金字塔。
CV_LKFLOW_PYR_B_READY 在调用和存储到 pyrB之前，先计算第二帧
的金字塔。
CV_LKFLOW_INITIAL_GUESSES 在函数调用之前，数组 B已包含特征点
的初始坐标值。
【331~332】
这些标志特别是在处理视频时会用到。图像金字塔的计算量较大，所以应该尽可能
避免重复计算金字塔。计算得到的图像对的后面一帧被作为下次计算的图像对的初
始帧。如果调用者为函数分配了缓存（而不是让函数在内部分配），则函数返回时，
每幅图像的金字塔将被存储在其中。告诉函数金字塔已经被计算，则函数不会再进
行重复计算。类似地，如果从前一帧中计算得到点的运动，则在计算这些点在下一
帧中的位置就有了一个好的初始估计。
应用的基本过程很简单：输入图像，在 featuresA 中列出需要跟踪的点，然后调
用函数。函数返回后，检查status
数组以确定哪些点被成功跟踪，再检查
featuresB得到这些点新的位置。
现在我们再来讨论前面未讨论的问题：如何确定哪些特征点是易于被跟踪的。前面
介绍过 OpenCV 函数 cvGoodFeaturesToTrack(),这个函数使用 Shi 和 Tomasi
提出的方法可靠地解决了这个问题。在多数情况下，使用 cvGoodFeatures-
ToTrack()和 cvcalcopticalFlowPyrLK()的组合可以获得很好的结果。当然，读
者可以使用自己的准则来确定被跟踪的点。
下面来看一个简单的使用cvGoodFeaturesToTrack()和 cvCalcopticalFlow-
PyrLK()的组合的例子（例10-1);也可参考图10-10.
Flow vectors
图10-10:金字塔 Lucas-Kanade稀疏光流。中间的图像是左边图像的下一帧
视频；右边的图像表示计算出的“易于跟踪特征点”的运动（右下图为了增加可
视性在黑色背景中标识出光流矢量）
例10-1:金字塔 Lucas-Kanade 光流代码
／／ Pyramid L-K optical flow example
／／
＃include <cv.h>
＃include <cxcore.h>
＃include <highgui. h>
const int MAX_CORNERS =500;
int main(int argc, char** argv) {
／／ Initialize, load two images from the file system, and
／／ allocate the images and other structures we will need for
／／ results.
／／
IplImage* imgA=cvLoadImage ("image0.jpg",CV_LOAD_IMAGE_GRAYSCALE) ;
IplImage* imgB=cvLoadImage ("image1. jpg",CV_LOAD_IMAGE_GRAYSCALE);
CvSize img_sz = cvGetSize( imgA );
int win_size=10;
IplImage* imgC = cvLoadImage(
＂。。／Data/OpticalFlow1. jpg",
CV_LOAD_IMAGE_UNCHANGED
）；
／／ The first thing we need to do is get the features
／／ we want to track.
／／
IplImage* eig_image = cvCreateImage ( img_sz,IPL_DEPTH_32F,1);
IplImage* tmp_image = cvCreateImage( img_sz, IPL_DEPTH_32F,1 );
int corner_count =MAX_CORNERS;
CvPoint2D32f* cornersA = new CvPoint2D32f[MAX_CORNERS ];
cvGoodFeaturesToTrack(
imgA,
eig_image,
tmp_i.mage,
cornersA,
＆corner_count,
0.01,
5.0,
0,
3,
0,
0.04
）；
cvFindCornerSubPix(
imgA,
cornersA,
corner_count,
cvSize(win_size,win_size),
cvSize(-1,-1),
cvTermCriteria (CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,20,0.03)
）；
／／ Call Lucas Kanade algorithm
／／
char features_found[ MAX_CORNERS ];
float feature_errors[ MAX_CORNERS ];
CvSize pyr_sz = cvSize( imgA->width+8, imgB->height/3 );
IplImage* pyrA = cvCreateImage( pyr_sz, IPL_DEPTH_32F, 1 );
IplImage* pyrB = cvCreateImage( pyr_sz,IPL_DEPTH_32F,1);
CvPoint2D32f* cornersB = new CvPoint2D32f[ MAX_CORNERS ];
cvcalcoptica1FlowPyrLK(
imgA,
imgB,
pyrA,
pyrB,
cornersA,
CornersB,
corner_count,
cvSize( win_size,win_size ),
5,
features_found,
feature_errors,
cvTermCriteria( CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20,.3 ),
0
）；
／／ Now make some image of what we are looking at:
／／
for( int i=0; i<corner_count; i++){
if( features_found[i]==0|| feature_errors[i]>550 ){
printf("Error is &f/n", feature_errors [i]);
continue;
｛
printf("Got it/n");
CvPoint po = cvPoint(
cvRound( cornersA[i].x ),
cvRound( cornersA[i].Y )
）；
CvPoint p1 = cvPoint(
cvRound( cornersB[i].x),
CVRound( cornersB[i].Y)
）；
cvLine( imgC, p0, p1, CV_RGB(255,0,0),2 );
｛
cvNamedWindow("ImageA",0);
cvNamedWindow("ImageB",0);
cvNamedWindow("LKpyr_Opt icalFlow" , 0) ;
cvShowImage("ImageA" , imgA) ;
cvShowImage ( "ImageB" , imgB);
cvShowImage("LKpyr_OpticalFlow",imgC);
cvWaitKey(0);
return 0;
稠密跟踪方法
OpenCV 中还有另外两种光流方法，但现在已经很少用到了。这些函数比 Lucas-
Kanade 方法慢很多；另外，它们不支持图像金字塔匹配而不能用于跟踪大幅度的
运动。本节我们简单讨论一下这些函数。
【332~334】
Horn-Schunck方法
Horn 和 Schunck于1981年提出这种方法。此方法是首次使用亮度恒定假设和推导
出基本的亮度恒定方程的方法之一。Horn 和 Schunck 求解方程方法是假定一个速
度v,和v,的平滑约束。该约束是通过对光流速度分量的二阶导数进行规则化获
得：
0=(7+x7+x)27
0=(7+ay+a7)77-e@
这里α是不变的权重系数，称为规则化常数（regularization constant).α的值较大可
）
以获得更平滑（也就是更局部一致）的运动流向量。这是强制进行平滑的一个简单的
约束，其效果是惩罚光流变化剧烈的区域。与 Lucas-Kanade 算法一样，Horn-
Schunck 方法也要通过迭代来解微分方程。
void cvCalcopticalFloWHS(
const CvArr*
imgA,
const CvArr*
imgB,
int
usePrevious,
CvArr*
velx,
CvArr*
vely,
double
lambda,
CvTermCriteria
criteria
）；
【335~336】
其中 imgA和imgB必须是8位的单通道图像。velx和vely存放x和y方向的速
度，它们必须是32位的浮点数单通道图像。参数 usePrevious 指定算法使用从前
一帧计算的 velx和 vely 速度作为计算新的速度的初始值。参数 lambda 是与
Lagrange 乘子相关的权重。读者可能会问：“什么是Lagrange 乘子？”Lagrange
乘子出现在当我们（同时）最小化运动－亮度方程和平滑方程时；它表示当最小化时
赋予每一个方程误差的相对权重。
块匹配方法
读者可能会想：“光流到底有什么用？只是匹配上一帧和下一帧的像素而已。”有
一些人的确是这样用的。“块匹配”是对一类将图像分割成被称为“块”
［Huang95;Beauchemin95]的小区域的相似算法的统称。典型的块是正方形的，包
含了一定数目的像素。这些块可以重叠，实际中它们通常是重叠的。块匹配算法将
前一帧图像和当前图像划分成小块，然后计算这些块的运动。这一类算法在视频压
缩算法和计算机视觉的光流中扮演重要角色。
块匹配算法对像素的集合进行处理而非单个像素，所以返回的“速度图像”通常比
输入图像的分辨率低，但也有例外，这取决于块之间的重叠程度。下面的公式给出
了结果图像的尺寸：
①
这种情况下最好忽略这里的描述而设置1ambda为1.
W result=
Weo-Woec+
H resulr
Hoo-Haoe+Haee
ftoor
OpenCV 中的实现算法用螺旋搜索，找出（前一帧图像的）原始块的位置，然后再将
候选的新块同原始块进行比较。比较的结果是像素绝对差值的和（也就是L1距
离）找到很好的匹配后就终止搜索。函数原型如下：
【336】
id cvcalcopticalFlowBM(
const CvArr*
prev,
const CvArr*
curr,
CvSize
block_size,
CvSize
shift_size,
CvSize
max_range,
int
use_previous,
CvArr*
velx,
CvArr*
vely
）；
参数表示的意思简明易懂。参数prev和 curr 是前一帧图像和当前图像，它们都
是8位单通道的图像。block_size 是块的尺寸，shift_size 是块间移动的步长
（这个参数决定了块与块之间是否重叠以及如果重叠，重叠的程度又如何）。参数
max_range
是一个块周围邻域的尺寸，函数在下一帧中搜索此邻域以找到对应的
块。如果设置了参数use_previous,则 velx 和 vely 的值将作为块搜索的初始
值?.参数 velx和 vely是存储计算得的块的运动的32位的单通道图像。如前文
所说，运动的计算是在块级上进行的，所以输出结果的图像的坐标是对块而言的
（就是说，像素的集合体），而不是原始图像中的每个像素。
mean-shift和
mshift 跟踪
这一节中将要介绍两种跟踪方法，mean-shift 和 camshift(“camshaift”是
“continuously adaptive mean-shift”的缩写）。前者是可用于多种应用的通用的数据
①
如果use_previous == 0,块搜索区域与原始块的位置的距离为 max_range.如果
use_previous !=0,则块搜索区域的中心将被偏移Δx=vel(x,y)和
Δy=vel,(x,y).
分析方法（在第9章的分割问题中讨论过），计算机视觉正是这些应用之一。在介绍
mean-shift 的基本原理之后，会介绍 OpenCV 是如何用它在图像中进行跟踪运动
的。后一种方法camshift是建立在mean-shift之上，它可以跟踪视频中尺寸可能产
生变化的目标。
mean-shift
mean-shift 算法①是一种在一组数据的密度分布中寻找局部极值的稳定的方法。若
分布是连续的，处理过程就比较容易，这种情况下本质上只需要对数据的密度直方
图应用爬山算法即可?.然而，对于离散的数据集，这个问题在某种程度上是比较
麻烦的。
【337】
这里用“稳定”这个形容词是指统计意义上的稳定，因为 mean-shift 忽略了数据中
的 outliers,即忽略远离数据峰值的点。mean-shift仅对数据局部窗口中的点进行处
理，处理完成后再移动窗口。
mean-shift算法的步骤如下。
（1)选择搜索窗口。
窗口的初始位置；
?窗口的类型（均匀、多项式、指数或者高斯类型）；
?窗口的形状（对称的或歪斜的，可能旋转的，圆形或矩形）；
?窗口的大小（超出窗口大小则被截去）。
（2)计算窗口（可能是带权重的）的重心。
（3)将窗口的中心设置在计算出的重心处。
（4)返回第（2)步，直到窗口的位置不再变化（通常会）?.
mean-shift 算法更加正式的描述是，mean-shift 算法与核密度估计的规则有关。
①
mean-shift是一个相当深奥的话题，我们在这里讨论的目的主要是使读者有个直观的印
象。参考Fukunaga[Fukunaga90]、Comaniciu 和 Meer[Comaniciu99]可以了解 mean-shift
的由来。
②
在这里使用“本质上”这个词是由于 mean-shift 会受尺度变化的影响。准确说就是：
mean-shift等价于先对连续分布用。mean-shift核进行卷积，然后再应用爬山算法。
③
迭代过程通常由最大迭代次数或者两次迭代中心变化的程度进行限制；虽然如此，迭
代过程最后都会收敛。
“核”是一个局部函数（例如高斯分布）。如果在足够的点处有足够合适的带权重和
尺度的核，数据的分布便可以完全根据这些核来表示。与核密度估计不同，mean-
shift 只估计数据分布的梯度（变化的方向）。若变化为0的地方则表示是这个分布的
峰值（虽然可能是局部的）。当然在附近或其他尺度上可能还有峰值。
图10-11是mean-shift算法中的计算公式。
xx
Start with akernel K(x-X)=ck
approximation of a probabllty
dsrulbuton P()=1/1
K(x-x).Focusor,the greadient VP(x)=/≥K(x-x).
Let:g(x)=-k'(x),the derivative of the kernel and we get:
VP(2)=K1/8
Meanshift
vector
图10-11:mean-shift等式及其含义
这些计算公式可用一个矩形的核进行简化，将 mean-shift 矢量等式简化为计算图
像像素分布的重心：
x.=M
y=Ma
这里，零阶矩的计算如下：
＝00W
ΣΣI(x,y)
一阶矩的计算为：
【338】
M0=ΣΣx(x,2)
和
Mo1=ΣΣyI(x,2)
①
矩形核并不随着到中心的距离下降，而是一个突然变成零的突然转换。这个与高斯核
的指数衰减不同，与Epanechnikov核的随着到中心的距离的开方衰减也不同。
mean-shift 矢量告诉我们如何将 mean-shift 窗口的中心重新移动到由计算得出的此
窗口的重心的位置。很显然，窗口的移动造成了窗口内容的改变，于是我们又重复
刚才重新定位窗口中心的步骤。窗口中心重定位的过程通常会收敛到mean-shift矢
量为0(也就是，窗口不能再移动）。收敛的位置在窗口中像素分布的局部最大值（峰
值）处。由于“峰值”本身是一个对尺度变化敏感的量，所以窗口大小不同，峰值
的位置也不一样。
图 10-12 是一个二维数据分布及初始化窗口（此例中，窗口为矩形）的例子。箭头表
示迭代过程最终收敛到分布的局部峰值。正如前文所说的，我们可以看到由于处于
mean-shift 窗口外部点不影响 mean-shift的收敛性，所以寻找峰值的这一过程在统
计意义上是稳定的。
图10-12:mean-shift算法的过程。将初始窗口放置在数据点的二维数组上，
连续重定位（窗口）中心到数据分布的局部峰值处，直到收敛
1998年，人们认识到这种寻优的算法可以用于跟踪视频中的运动物体［Bradski98a;
Bradski98b],此后这个算法被大大地扩展［Comaniciu03].OpenCV 中运用 mean-
shift 算法的函数是在图像分析的层面上实现的。这也就是说，OpenCV 中实现的
mean-shift 把一幅代表将要被分析的密度分布的图像作为输入，而不是一组任意的
点（可能是任意的维数）作为输入。
【339~340】
int cvMeanShift(
const CvArr*
prob_image,
存在mhi中的最大持续时间。换句话说，mhi中任何比timestamp 减去 duration
的值早（少）的像素将被置为0.
一旦运动模板记录了不同时间的物体轮廓，就可以用计算 mhi 图像的梯度来获取
全局运动信息。计算出的这些梯度（例如，用第6章讨论过的 Scharr 或 Sobel 梯度
函数来计算），一些梯度值会很大并且是无效的。mhi 图像中被置0的旧的或没有
运动的部分的梯度是无效的，因为它们会在轮廓的外部边缘区域人为地产生很大的
梯度值，如图10-15(A)所示。由于使用cvUpdateMotionHistory()将新的轮廓引
进 mhi的时间步长是已知的，所以梯度值应该为多大也就知道了（即为 dx和 dy
step derivatives).那么我们可以用梯度幅值来消除特别大的梯度，如图10-15(B)所
示。最终，我们就得到全局运动的测量，如图 10-15(C).图 10-15(A)和（B)是由
cvCalcMotionGradient()函数得到的：
【343】
void cvCalcMotionGradient(
const CvArr*
mhi,
CvArr*
mask,
CvArr*
orientation,
double
deltal,
double
delta2,
int
aperture_size=3
）；
（B)
C
图10-15:mhi图像的运动梯度：（A)梯度值和方向；（B)大的梯度被去除；（C)
得到全局运动方向
函数cvCalcMotionGradient 中，所有的图像矩阵都是单通道的。函数的输入
mhi是一幅浮点数值的运动历史图像，delta1和delta2则分别是允许的最小和最
大梯度值。期望的梯度值是连续调用 cvUpdateMotionHistory()的每个轮廓间当
前时间标记的平均数；将delta1 设置为小于此平均值的一半，delta2 设置为大
于此平均值的一半比较合适。变量 aperture_size 设置梯度算子的宽和高，其值
可以为－1(3x3 CV_SCHARR梯度滤波器）、3(默认的3x3 Sobel滤波器）、5(5x5 Sobel
滤波器）或者7(7x7滤波器）。函数的输出mask是一幅8位的单通道图像，其中的
378 第10章
CvRect
window,
CvTermCriteria
criteria,
CvConnectedComp*
comp
）；
CVMeanShift()函数中，参数 prob_image是单通道图像，代表可能位置的密度，
其类型可以是 byte 或 float.参数 window的位置设置在为指定初始位置，大小为
核窗口的大小。终止条件criteria 主要由 mean-shift移动的最大迭代次数和可视
为窗口位置收敛的最小移动距离组成，其定义在其他章节中有描述①.连接部件
comp 在comp->rect中包含收敛后搜索窗口的位置，在comp->area中包含了窗口
内部所有像素点的和。
函数cvMeanShift()使用的是矩形窗口的mean-shift 算法，但这同样可以用于跟
踪。在这种情况下，需要首先选择代表物体的特征的分布（例如，颜色＋纹理），然
后在物体的特征分布上开始mean-shift窗口搜索，最后计算下一帧视频中所选择的
特征的分布。从当前的窗口位置开始，mean-shift 算法寻找特征分布的新的峰值或
者 mode,它们（被假定）设置在最初产生颜色和纹理的物体的中心。这样，mean-
shift窗口就可以逐帧跟踪物体的运动。
CamShift
另一个相关的算法是CamShift跟踪器。与mean-shift 不同的是，CamShift 搜索窗
口会自我调整尺寸。如果有一个易于分割的分布（例如保持紧密的人脸特征），此算
法可以根据人在走近或远离摄像机时脸的尺寸而自动调整窗口的尺寸。CamShift
算法的形式如下：
?int cvCamShift(
const CvArr*
prob_image,
CVRect
window,
CvTermCriteria
criteria,
CvConnectedComp* comp,
CVBox2D*
box = NULL
）；
前四个参数的含义与cVMeanShift()算法中的含义相同。如果参数box 不为空，
①
强调一下mean-shift 总是会收敛，但是如果一个分布在局部峰值附近非常“平坦”，
收敛将会变得非常慢。
Current
TOIE
SIllouette
Current
Current
Sllouette
Sillouette
Current
则它包含了新的尺寸的box,其中也包含根据二阶矩计算出的物体的方向。在跟踪
的应用中，会把由前一帧计算出的新尺寸的box作为下一帧的window.
注意：很多人认为mean-shift 和 camshift 利用颜色特征进行跟踪，但实际上不完全是这
样。两种算法跟踪在prob_image 中所代表的任意的分布，所以它们是轻量级的、稳
定和有效的跟踪器。
运动模板
运动模板（motion template)是由 MIT 媒体实验室的 Bobick 和 Davis 提出的
［Bobick96;Davis97],并且其中的一人又参与了它进一步的开发［Davis99;
Bradski00].而这一步工作为运动模板在OpenCV 中的实现提供了基础。运动模板
是一种有效的跟踪普通运动的方法，尤其可应用在姿态识别中。运用运动模板需要
知道物体的轮廓（或者轮廓的一部分，此处轮廓指 silhouette,即实心轮廓），而轮廓
的获取有不同的方法。
【341~342】
（1)最简单的获取物体轮廓的方法是用一个静止的摄像机，再使用帧间差（见第9
章论述）得到物体的运动边缘，这种信息就足够让运动模板发挥效用。
（2).利用颜色信息。例如，如果已知背景的颜色，比如亮绿色，那么可以很简单地
将不是亮绿色的任何物体看作前景。
（3)另一种方法（在第9章中讨论过）是学习一个背景模型，从此背景中可以将新的
前景物体／人以轮廓的形式分割出来。
（4)使用主动轮廓技术。例如，创建一个近红外光的墙，将能感应近红外线的摄像
机对着墙，那么任何介于墙与摄像机之间的物体都会以轮廓的形式出现。
（5)使用热感图像。任何温度高的物体（如人脸）都可以被认为是前景。
（6)最后，可以用第9章中谈到的分割技术（例如金字塔分割或 mean-shift分割）来
生成轮廓。
从现在开始，假设我们有一个如图10-13(A)中白色矩形所代表的被较好的分割出的
物体轮廓。白色表示所有像素都被设置为最新的系统时间的浮点数值。随着矩形的
运动，新的轮廓将被捕获并且被（新的）当前轮廓覆盖；新的轮廓为图 10-13(B)和
图 10-13(C)中的白色矩形所示。较早的运动在图10-13中表示为亮度渐暗的矩形。
这些连续变暗的轮廓记录了早前运动的历史，所以被称为“运动历史图像（motion
history image)”。
图10-13:运动模板图：（A)当前时间分割出的物体（白色）；（B)在下一个时间
点，物体运动并以（新的）当前时间标记，将前面的分割边界留在后面；（C)在下
一个时间点，物体继续运动并将此前的分割标记为渐暗的矩形，这些包含运动
的序列就产生运动历史图像
【342】
若轮廓的时间超过设定的比当前时间早的持续时间，则其被置0,如图10-14所
示。OpenCV中完成运动模板构建的函数是cvUpdateMotionHistory():
void cvupdateMotionHistory(
const CvArr*
silhouette,
CvArr*
mhi,
double
timestamp,
double
duration
）；
（A)
（B)
运动
运动
图10-14:两个物体的运动模板轮廓（A);超过设定持续时间的轮廓被置0(B)
在cvupdateMotionHistory()中，所有图像都是单通道图像。图像
silhouette
是一幅单字节图像，其中非0像素代表前景物体的最新的分割轮廓。图像 mhi是
一幅浮点值的图像，代表运动模板（也就是运动历史图像）。timestamp 是当前系统
时间（典型地以毫秒为单位），duration 如刚才谈到的表示运动历史像素的允许保
存在mhi中的最大持续时间。换句话说，mhi中任何比timestamp 减去 duration
的值早（少）的像素将被置为0.
一旦运动模板记录了不同时间的物体轮廓，就可以用计算 mhi 图像的梯度来获取
全局运动信息。计算出的这些梯度（例如，用第6章讨论过的 Scharr 或 Sobel 梯度
函数来计算），一些梯度值会很大并且是无效的。mhi 图像中被置0的旧的或没有
运动的部分的梯度是无效的，因为它们会在轮廓的外部边缘区域人为地产生很大的
梯度值，如图10-15(A)所示。由于使用cvUpdateMotionHistory()将新的轮廓引
进 mhi的时间步长是已知的，所以梯度值应该为多大也就知道了（即为 dx和 dy
step derivatives).那么我们可以用梯度幅值来消除特别大的梯度，如图10-15(B)所
示。最终，我们就得到全局运动的测量，如图 10-15(C).图 10-15(A)和（B)是由
cvCalcMotionGradient()函数得到的：
【343】
void cvCalcMotionGradient(
const CvArr*
mhi,
CvArr*
mask,
CvArr*
orientation,
double
deltal,
double
delta2,
int
aperture_size=3
）；
（B)
C
图10-15:mhi图像的运动梯度：（A)梯度值和方向；（B)大的梯度被去除；（C)
得到全局运动方向
函数cvCalcMotionGradient 中，所有的图像矩阵都是单通道的。函数的输入
mhi是一幅浮点数值的运动历史图像，delta1和delta2则分别是允许的最小和最
大梯度值。期望的梯度值是连续调用 cvUpdateMotionHistory()的每个轮廓间当
前时间标记的平均数；将delta1 设置为小于此平均值的一半，delta2 设置为大
于此平均值的一半比较合适。变量 aperture_size 设置梯度算子的宽和高，其值
可以为－1(3x3 CV_SCHARR梯度滤波器）、3(默认的3x3 Sobel滤波器）、5(5x5 Sobel
滤波器）或者7(7x7滤波器）。函数的输出mask是一幅8位的单通道图像，其中的
非0值代表此点处的梯度有效，orientation 是一幅浮点值的图像，给出每一点
梯度方向的角度。
函数cvCalcGlobalorientation()计算有效梯度方向矢量和来获得全局运动
方向。
double cvCalcGlobalOrientation (
const CvArr*
orientation,
const CvArr*
mask,
const CvArr*
mhi,
double
timestamp,
double
duration
）；
函数的输入为由cvCalcMotionGradient()计算得到的 orientation和mask图像
以及timestamp、duration和 cvupdateMotionHistory()得出的mhi;函数的输
出为如图 10-15(c)所示的矢量和的全局方向。timestamp和 duration给函数设定
了mhi和 orientation 图像中需要考虑的运动量。用每个mhi轮廓的重心可以计
算得到全局运动，但是对已计算出的运动矢量求和会快很多。
我们还可以将运动模板 mhi 图像分割为独立的区域，并计算区域里的局部运动，
?如图10-16所示。图中，搜索mhi图像寻找当前的轮廓区域。当找到被标记为最新
当前时间的区域时，继续搜索区域的边界以寻找紧邻边界外围的最新近的运动（最
近的轮廓）。若找到这样的运动，就会用逐级洪水泛滥法将从感兴趣物体的当前位
置分割出来，计算分割区域局部运动的梯度方向，再移除此区域。重复此过程直到
所有的区域都被找到（如图10-16所示）。
【344~345】
分割和计算局部运动的函数为cvSegmentMotion():
CvSeg* cvSegmentMotion(
const CvArr*
mhi,
CvArr*
seg_mask,
CvMemStorage*
storage,
double
timestamp,
double
seg_thresh
）；
函数中，mhi 是单通道、浮点值的输入图像。我们还传入 storage,它是一个由
cvCreateMemStorage()分配的 CvMemoryStorage 结构。另一个输入为
timestamp,它的值是要从其中分割出局部运动的 mhi 图像中最新轮廓的值。最
后，还必须传入一个参数seg_thresh,它是认为可以接受为相关联运动的的最大
从当前时间到早先的运动。之所以要输入这个参数，是由于最近的和早前的运动的
轮廓可能存在重叠，但你并不想将它们连接在一起。
（A)
（B)
（C)
（D)
图10-16:分割mhi图像中的局部运动区域：（A)搜索mhi图像寻找当前轮廓
a.当搜索到此轮廓，沿其边界搜索其他较近的轮廓（b);若找到，用逐级洪水
泛法c分割出局部运动；（B)计算分割出的局部运动区域的梯度，并用此梯度
计算局部运动；（C)移除此区域并搜寻下一个当前轮廓的区域d,沿着其搜索
（e),再用逐级洪水泛滥法填充f;(D)计算新分割区域内的运动，重复（A)(C)过
程直到没有剩余的当前轮廓
【345】
通常最好将 seg_thresh的值设置为轮廓时间平均差别的1.5倍。函数返回
CvConnectedComp 结构的 CvSeq 序列，其中每一个代表搜索到的一个描述局部运
动区域的分割的运动；函数还返回一个单通道的浮点值图像，seg_mask,其中每
个被分割出的运动的区域被标记为不同的非0值（seg_mask中值为0的像素表示无
运动）。用从CvConnectedComp或者 seg_mask 中特定值中选择的掩码区域，调用
cvCalcGlobalorientation()一次来计算一个局部运动，如下所示：
cvCmps(
seg_mask,
／／［value_wanted_in_seg_mask] ,
／／［your_destination_mask],
CV_CMP_EQ
讨论到这里，现在应该能够理解OpenCV的。。。／opencv/samples/c/目录中的
motempl.c例子。我们现在从 motempl.c中 update_mhi()函数抽出一些关键点来解
释。update_mhi()函数使用有阈值的帧间差来获取模板，然后得到的轮廓传递给
cvupdateMotionHistory():
?..
cvAbsDiff( buf[idx1], buf[idx2], silh );
CvThreshold( silh, silh, diff_threshold, 1, CV_THRESH_BINARY );
cvUpdateMotionHistory( silh, mhi, timestamp, MHI_DURATION );
。。。
计算得到的 mhi 图像的梯度，再由 cvCalcMotionGradient()获得有效梯度的掩
码。分配CvMemStorage(或者，如果已经存在则释放），分割出的局部运动存储在
cvConnectedComp结构中，而cvConnectedComp在包含CvSeq结构的seq中：
?.
CvCalcMotionGradient(
mhi,
mask;
orient,
MAX_TIME_DELTA,
MIN_TIME_DELTA,
3
）；
if(!storage)
storage =cvCreateMemStorage(0);
else
cvClearMemStorage(storage) ;
seq =cvSegmentMotion(
mhi,
segmask,
storage,
timestamp,
MAX_TIME_DELTA
）；
【346~347】
我们用一个for循环重复从sep->total个CvConnectedComp结构中提取每个运动
的边界矩形。循环从－1开始，这样指定为的是要寻找整幅图像的全局运动。对于
局部运动块，忽略面积小的分割后，用 cvCalcGlobalOrientation()计算方向。
函数将运动的计算限制在包含局部运动的感兴趣区域（ROI),而不是提取掩码；然
current
current
Jiouere!
Sillouette
Slllouette
Current o
（
后，计算在局部 ROI 中有效运动的实际位置。任何太小的这样的运动区域将被丢
弃。最后，函数将运动绘制出来。图10-17表示了对一个挥动手臂的人使用此函数
的输出结果，即为两行每行四帧连续原始图像上方的图像。（完整的代码可以参
考．．．／opencv/samples/c/motempl.c.)序列中，“Y”姿态由第8章讨论的形状算子
（Hu矩）识别得出。samples代码中还未包含形状识别。
Y pose
Y 0050
图10-17:运动模板函数运行结果。从左至右、从上到下，一个人在运动，输
出结果中大的八边形表示全局运动，小的八边形表示局部运动；“Y”姿态可
由形状算子（Hu矩）识别
for( i =-1; i <seq->total; i++){
if( i<0){// case of the whole image
／／。。。［does the whole image]...
else{//i-th motion component
comp_rect =((CvConnectedComp*)cvGetSeqElem( seq, i ))->rect;
／／ [reject very small components]...
（
。。。［set component ROI regions]...
angle = cvcalcGlobalOrientation( orient, mask, mhi,
timestamp,MHI_DURATION);
。。。［find regions of valid motion]. . .
。。。［reset ROI regions]...
。。。［skip small valid motion regions]. . .
。。。［draw the motions]...
【347】
预估器
假设我们跟踪一个穿过摄像机视场的人，每一帧都要确定此人的位置。正如我们知
道的，有多种方法可以达到这个目的，但是我们发现每种情况下得到的只是此人在
每帧中位置的估计。这个估计，不太可能相当的准确。造成不准确的原因有很多。
可能是传感器的不精确，早期处理阶段的近似，遮挡或者阴影，又或者是由于腿和
手臂在行走中的摆动造成的明显的形状变化。不管源头是什么，我们期望这些测量
对可能从理想传感器获得的“实际”值会发生变化，而且可能是任意的。整体考虑
这些影响因素，可以将这些不准确简单地看作向跟踪过程添加噪声。【348~349】
我们希望可以最大限度地使用测量结果来估计此人的运动。所以，多个测量的累积
可以让我们检测出不受噪声影响的部分观测轨迹。一个关键的附加要素即此人运动
的模型。例如，我们可能用下面的方式建立此人的运动模型：“一个人从图像的一
边进入并以恒定的速度穿过图像。”有了这个模型，我们不仅可以知道此人在什么
位置，同时还可以知道我们的观察支持模型的什么参数。
这个任务分成两个阶段（如图10-18 所示）。在第一阶段，即预测阶段，用从过去得
到的信息进一步修正模型以取得人（或物体）的下一个将会出现的位置。在第二阶
段，即校正阶段，我们获得一个测量，然后与基于前一次测量的预期值（即模型）进
行调整。
完成两个阶段估计任务的方法被归入预估器（estimator)这个标题下，其中 Kalman
滤波器［Kalman60]是最广泛使用的。除了 Kalman 滤波器，另一重要的方法是
condensation 算法，它是更为广泛的一类方法粒子滤波在计算机视觉领域的实现。
Kalman 滤波器和 condensation 算法主要的区别在于状态概率密度是如何描述的。
后面将探究这个区别的含义。
【349~350】
预测阶段“根据前面数据预测
校正阶段“根据新测量调整
图 10-18:两阶段的预估器循环：基于前面数据的预测和根据最新测量的调整
Kalman 滤波器
Kalman 滤波器（即卡尔曼滤波器）最初于1960年提出，之后发展成为信号处理领域
中不同方面的重要方法。Kalman 滤波器的基本思想是，若有一组强而合理的假
设，给出系统的历史测量值，则可以建立最大化这些早前测量值的后验概率的系统
状态模型?.Welsh 和 Bishop[Welsh95]做了比较好的介绍。另外，无需存储很长的
早前测量历史，我们也可以最大化后验概率，即重复更新系统状态模型，并只为下
一次更新保存模型。这样就大大地简化了这个方法的计算机实现。
在详细介绍这些在实际应用中的含义之前，我们花一点时间来看看提到的假设是什
么。理论上，Kalman 滤波器需要三个重要的假设：（1)被建模的系统是线性的；（2)
影响测量的噪声属于白噪声；（3)噪声本质上是高斯分布的。第一条假设的意思是k
时刻的系统状态可以用某个矩阵与k-1时刻的系统状态的乘积表示。余下两条假
①
“合理的”在这里的意思是“限制非常宽松使得这个方法对真实世界中出现的相当多
的实际问题都有用”。“合理的”只是仅次于“完美的”。
②
修饰词“后验的”是学术领域用于描述“事后解释”的一个术语。当我们说这样一个
分布“最大化后验概率”时，意思是虽然本质上是对“实际发生过的”一种可能的解
释，但结合已经观测到的数据，这个分布实际上是最可能的一个，也就是以回顾的方
式分析。
3
设，即假设噪声是高斯分布的白噪声，其含义为噪声与时间不相关，且只用均值和
协方差（也就是噪声完全由一阶矩和二阶矩描述）就可以准确地为幅值建模。这些假
设看起来似乎非常苛刻，但它们实际上却应用在非常广泛的普通环境中?.
“最大化前期测量值的后验概率”是什么意思？意思就是在获得测量值－考虑早
前的模型和新测量值的不确定性－之后新建立的模型是具有最高正确概率的模
型。给定三条假设，Kalman 滤波器是将从不同来源获得的数据或从同一来源不同
时间获得的数据结合的最好的方法。从我们知道的信息开始，获取新的信息，然后
根据对旧信息和新信息的确定程度，用新旧信息带权重的结合对我们知道的信息进
行更新。
【350~351】
我们需要用一点数学和一维运动的情况来说明上文所讲的概念，你可以跳过下一
节。线性系统和高斯分布是非常易懂的，如果你看都不看一眼，Kalman 博士可能
会不开心。
Kalman滤波器相关的一些数学知识
Kalman 滤波器的核心是什么？是信息融合。假设你想知道某一点在一条直线的什
么位置（这就是一维的情况）?,由于噪声的影响，将会得到两个不可靠的（从高斯概
率的角度）物体位置的变量：x1和x2.因为这些测量值中存在高斯不确定性，所以
有均值x和x以及标准差σ,和σ2.标准差其实是反映了我们对于测量值好坏程度
的不确定性。以位置为变量的概率分布为下面的高斯分布：
a./2x(-(2)(=1.2)
p,(x)=
若有这样两个具有各自高斯概率分布的测量，我们期望在同时给出两个测量的条件
下某一x值的概率密度与p(x)=p1(x)p2(x)成比例。这个乘积的结果也是一个高斯分
布，可以按照下面的方式来计算这个新的分布的均值和标准差。由于有
①
再多说一句。在这里我们实际上又用了另一个假设，即初始分布在本质上也是高斯分
布的。通常应用中的初始状态是准确地知道的，或者至少可以这样处理，于是这就满
足了我们的要求。如果初始状态为（举例说）在卧室或者在浴室的概率各占一半，则我
们会很不幸运，则需要比单一的Kalman滤波器更有好的方法。
②
更多详细的解释以及一个曲线轨迹的例子，可以参考 J.D.Schutter,J.De Geeter,
T.Lefebvre 和 H.Bruyninckx
所著的“
Kalman Filters: A Tutorial”
（http://citeseer.ist.psu.edu/443226.html).
并且高斯分布在平均值处最大，我们可以简单地用计算关于x的p(x)的导数来获得
平均值。由于函数在其导数为0处值最大，所以有
由于概率分布函数p(x)不会为0,所以括号中的项必须为0.解这个等式后，可得
到一个非常重要的关系式：
（0+)+(+)=x
【351】
于是，新分布的均值即为两个测量到的均值的加权组合，而权重则由两个测量相关
的不确定性决定。可以看到，例如，如果第二个测量的不确定性σ2特别大，则新
的均值与更确定的前期测量的均值x,本质上是相同的。
知道了新的均值x2,将其代入P12(x)的达式中进行整理，则不确定性o/2即为：
σ22
12=
20+20
计算到这里，你可能想知道这个说明了什么。实际上，当我们用新的均值和不确定
性进行新的测量时，可以将新的测量与已有的均值和不确定性结合来获得由更新的
均值和不确定性确定的新的状态。（这些过程现在有数值表达式，我们马上会讲
到）。
两个高斯测量相结合等价于一个高斯测量（其均值和不确定性是可计算的）这个性质
对我们来说是最重要的一个。这个性质的意思是，当有M个高斯测量时，可以先
合并前两个，再将前两个合并的结果与第三个进行合并，然后再将合并的结果与第
四个合并，以此类推。计算机视觉中的跟踪就是这样使用的；由一个测量推出下一
个再推出其下面的一个。
①
整理的过程有些繁杂。如果你想验证整个过程，这样做会容易很多：（1)根据x12和012
由高斯分布P12(x)的等式开始，（2)将等式中的与x12相关的用x和x2代替，与O12相关
的用σ,和O2代替，最后（3)验证替换后的结果是否能分离成最开始的高斯的乘积。
由于测量值（x,,σ,)是时间离散的，可以如下面所述的方法来计算（x,6)估计值的当
前状态。在时间点1,只有第一个测量x=x,和它的不确定性＝02.在最优估计
等式中将它们替换，则得到一个迭代等式：
重新改写一下等式，则可得到如下有用的形式：
x2=x1+
（－3)+
在我们思考这个公式如何使用之前，先来计算一下2的另一种表达形式。首先，
令？？＝2,则有：
【352】
100
62=
0+29
像处理x,那样重新整理一下等式，则得到新的测量的估计协方差的迭代公式：
02=1+
用这些等式当前的形式，我们可以很清晰地分辨出“旧”信息（在新的测量结果得
到之前已知的）和“新”信息（最新的测量结果）。在第二个时间点 2,新的信息
（x2-x)叫“变化（innovation)”。现在我们也可以得到最优迭代更新比例：
0
K=
2+62
这个比例叫更新率（update gain).用K的这个定义，我们可以得到下面的比较方便
的递归形式：
x2=x1+K(x2-x1)
32=(1-K)2
在Kalman 滤波器文献中，如果讨论的是普通的测量，则第二个时间点“2”通常
就表示为k,而第一个时间点则为表示为k-1.
动态系统
在简单的一维情况的例子中，我们认为目标在某一点x,并且在这一点进行连续的
测量。在这种情况中，我们并没有特别地去考虑到在测量间物体实际上可能在运动
的情况。对于这种情况，我们使用预测。在预测阶段，在新的测量添加进来之前，
我们利用已知去计算系统的期望值。
实际上，预测阶段发生在新的测量完成之后和新测量的值被添加进系统状态的估计
中之前。举例来说，我们在时间t测量汽车的位置，再在时间t+dt测量一次。若汽
车的速度为v,则我们不直接合并第二次测量的结果。首先，用在时间 t得到的模
型向前推，得到系统在时间t+dt,即新的测量值被合并之前的时刻的模型。这样，
在时间 t+dt 获得的新的测量值与由旧模型向前映射到时间t+dt的模型融合，而不
是直接与这个旧的模型融合。这就是图10-18中描绘的循环所表示的意思。在
Kalman滤波器应用中，我们将要考虑三种运动。
【353~354】
第一种运动是动态运动（dynamical motion).这种运动是我们期望的前次测量时系
统状态的直接结果。如果我们在位置x、时间t时测量速度为v的系统，则在时间
t+dt 我们可以预期此系统在位置x+vxdt,速度可能仍为此速度。
第二种形式的运动叫控制运动（control motion).这种运动是我们期望的，由于某种
已知的外部因素以某种原因施加于系统。正如运动名称所暗示的，控制运动最普遍
的一个例子是，当对我们施加了控制的系统估计其状态时，我们知道我们的控制会
使系统产生什么样的运动。这种情况对机器人系统来说尤其是这样。在机器人系统
中，这种控制就是系统告诉机器人（例如）加速或者向前。很明显，在这种情况下，
如果机器人在时间t时在位置x并且运动速度为v,那么由于我们命令机器人加
速，在时间 t+dt我们预期机器人运动的距离不是x+vxdt(这个值是系统没有外在控
制时的预期值），而是比这更远。
最后一类运动是随机运动（random motion).即便是在简单的一维情况的例子中，如
果观测的目标有因任一原因而产生运动的可能性，那么就需要在预测阶段包含这种
随机运动。这种随机运动的影响是简单地增加状态估计随时间的协方差。随机运动
包含未知的或不在我们控制中的任何运动。但是，在Kalman 滤波器中，随机运动
被假设成为高斯模型（也就是一种随机）或者至少是可以用高斯模型有效地来建模。
因此，我们在融合一个新的测量值前先进行“更新”以便将动态变化包含进仿真模
型中。在更新阶段首先将我们所已知的目标的运动信息应用到其先前的状态上，然
后再应用任何由我们施加的控制或从某个外部代理施加的已知的控制得到的任何附
加信息，最后将在上次观测完成后可能改变系统状态的符合定义的随机事件融合进
来。应用这些因素之后，我们便可以融合下一个新的测量值。
实际上，在系统“状态”比我们的仿真模型复杂得多的时候，动态运动尤其重要。
通常，当目标在运动时，系统的“状态”有多个组成成分，例如位置和速度。当然
在这种情况下，系统状态随我们认为其具有的速度而发展。在下一节中将讨论如何
处理具有多个状态成分的系统。我们还将会用一些稍微复杂的符号来处理这些问题
的新方面。
【354】
Kalman 方程
现在我们来推广简单模型中的运动等式。这个更普遍的讨论中的模型是目标状态的
线性函数F.比如说，这样的一个模型可能会考虑到先前运动一阶和二阶导数的结
合。我们也将会看到如何在模型中处理控制输入uk.最后，我们将得到一个更符
合实际的观测模型z,在这个模型中只需要测量几个模型的状态变量，并且测量值
与状态变量间没有直接联系。
首先，我们看一下前一节中提到的更新率 K如何影响状态估计。若新的测量值的
不确定性非常大，则此新的测量值对状态估计不起作用，等式还原为我们已知的在
时间t-1的结果。相反地，如果开始时原始测量的协方差非常大，那么我们进行新
的更准确的测量，然后“相信”这个值。若两个测量值具有相同的确定性（协方
差），则新的期望值必将在它们之间。以上的论述与我们的合理的期望是相符的。
图10-19表示不确定性是如何根据新的观测随时间变化。
N(zy )
实测
现在知道的
过去知道的
（x)N
图10-19:将先验知识N(Xk-1,0k-1)与观测N(Zk,Ok)结合；结果为新的估计
（x)N
更新，对不确定性敏感，可以推广到其他状态变量。最简单的例子就是在视频跟踪
的应用中，视频中的物体可以具有两个或三个方向的运动。一般来说，状态可能包
①
观察从xk到zk在符号上的变化。后者是文献中的标准，意思是zk不仅是（有时是不是）
对位置xk的测量，可能是对模型多个参数的更一般的测量。
含附加的元素，例如被跟踪的物体的速度。在这些任意的一般情况下，为了说明我
们所要讨论的问题，需要再引入一些符号。我们将时间k的状态描述推广为时间
k-1的状态的函数。
【355~356】
＊＋＇ng+'*x4='x
这里，xx现在是一个状态元素的n维向量，F是一个与xx相乘的nxn矩阵，其有时
也被称作传递矩阵。向量 uk是新添加的。它的作用是允许外部控制施加于系统，
由表示输入控制的c维向量组成。B是一个联系输入控制和状态改变的nxc矩阵。
变量 wk是一个关联直接影响系统状态的随机事件或外力的随机变量（通常称为过程
噪声）。假设wk的元素具有高斯分布N(0,Qn),nxn协方差矩阵Q(Q可以随时间变
化，但通常不这样做）。
聪明的读者，或者已经对Kalman 滤波器有一定了解的读者会注意到漏掉的另一个
重要假设，也就是在控制 uk和状态变化间存在线性关系（通过矩阵乘法）。在实际
应用中，这条假设常常首先被违反。
一般来说，测量值 zk有可能是但也有可能不是状态变量的直接测量xk.(例如，如
果想要知道一辆汽车行驶得多快，可以用雷达枪测量它的速度，也可以测量排气管
传出的声音，在前一种情况中，zk为xk加上测量噪声，而后一种情况，这种关系
并不这么直接）。我们将这种情况总结为测量下面等式给出的测量值 zk的m维向
量：
＊a+*x*H='z
这里，Hk是mxn矩阵，Vk是测量误差，也被假设为具有高斯分布N(0,Rk)和mxm
协方差矩阵Rk.
在我们还没有彻底被弄糊涂之前，我们来看一个具体的测量在停车场行驶的汽车的
实际例子。我们可以将汽车的状态用两个位置变量x和y,以及两个速度变量vk和
V)表示。这四个变量组成状态向量xx的元素。这就是说F的正确形式为：
x
10dt 0
y
F=
01
0 dt
【356】
＝
＇x
V
00
10
0001
①
这些项中的k表示它们可以随时间变化，但是并不需要这样。在实际的应用中，H和
R通常不随时间变化。
但是，当使用摄像机去测量汽车的状态时，可能只能测量到位置变量：
ZK
这就暗示了H结构类似于下面这样：
1
0
0 1
H=
00
0
0
在这个例子中，我们可能并不真正认为汽车的速度是不变的，所以要设置一个值
Qk来反映这个问题。我们对视频流用图像分析方法来测量汽车的位置，再根据对
测量精确程度的估计来选择Rk.
现在剩下所要做的就是将这些表达式嵌入推广形式的更新等式中。基本思想是相同
的，但是，首先要计算一下状态的先验估计xi.在文献中，比较常见的（虽然不是
普遍的）是用符号作为上标来表示“在新的测量之前”；在此我们也采用这个惯
例。这个先验估计如下：
x=Fx1-1+BUk-1+Wk
用Pr表示误差协方差，此协方差在时间k的先验估计由其在时间k-1的值获得：
Pr=FP,FT+Q41
这个等式就构成了预估器预测部分的基础，它告诉我们根据我们已经看到的“我们
可以期望什么”。由此，我们可以给出（并没有引出）所谓的 Kalman 更新率或者混
合比例，它告诉我们如何给定新信息相对已经知道的信息的权重：
K=PTH[(HPTHT+R)-1
虽然这个等式看起来有些勉强，但实际上它也不是很恐怖。若思考一下不同的简单
的情形，就能更容易地理解这个等式。在直接测量一个位置变量的一维的例子中，
Hk仅仅是只有1的1x1矩阵！所以，如果测量误差是241,那么Rk也是只含有1
值的1x1矩阵。类似的，Pk也正是这个协方差0/.于是，那个庞大的等式就简化
为如下：
K=
02
【357】
σ2+02
跟踪与运动
391
第1页
OCR文字识别结果:
注意，这正是我们期望的。在前一节中首次提到的更新率让我们在获得新的测量值
时，可以计算xk和Pk最优的更新值
P., = (-K, H,)P
这些方程又一次让人乍一看很头疼，但是在所讨论的简单的一维情况中，它们也不
是这么槽的。最优权重和更新率以与一维情形相同的方法获得，除了此次我们在解
之前设置对x的偏导数为0来最小化位置状态x的不确定性。可以先设置F=八1是
单位矩阵)，B=O和Q=O,来表示更简单的一维情况。在更一般的等式中作如下的
替换，同时也揭示出与一维滤波器导数的相似之处x,4一毛，x,呓一元，K,〈-K,
OpencV和Kalman滤波器
根据我们知道的这一切，你可能会觉得不需要OpencV为我们做什么事情或者迫
切地需要OpencV为我们完成一切工作。幸运的是，OpencV司'以满足两种需求。
OpencV提供了四个与应用Kalman滤波器直接相关的函数。
cvcreate Kalman(
int n Dynamparams,
int nMeasureparams,
int ncontro1Params
cvRelease Kalman(
CvKalman** kalman
第一个函数产生和返回一个cvKalman数据结构的指针，第二个函数删除这个
结构。
typedef struct CvKalman
int MP // measurement vector dimensions
int DP // state vector dimensions
int CP // control vector dimensions
CvMat* state_pre ; // predicted state
CvMat* state_post // corrected state
CvMat* transition_matrix // state transition matrix
/ / F
392第10章
第2页
OCR文字识别结果:
冒匿
CvMat* control_matrix // control matrix
f / B
// ( not used I. f there is no control )
CvMat · measurement_mat r ix
. / / measurement mat ri x
// H
CvMat* process_noise_cov
// process noise covariance
!l Q
cvMat · measuremnt-noise-cov
/ / measurement noise covariance
CvMat · error-cov-pre
CvMat* gain // Kalman
CvMat * error-covJost
CvMat * tempi
CvMat* temp2
CvMat* temp3
CvMat* temp4i
CvMat * temps
/ / R
/ / prior'error covariance
/ / ( P-k'. F P-k-1 Ft ) + Q
gain matrix
/K一大·P一大'H^T(H P一大'H八T+R)。一二
/ / posteriori error covariance
/P一大，(工一K一k H)P一大'
/ / temporary matrices
CvKalman
[358-359]
下面两个函数是Kalman滤波器的实现。一旦cvKalman结构被赋值，则调用
cvKalmanpredict O计算下一个时间点的预期值，然后调用cvKalmancorrect 0
校正新的测量值。这些函数每个运行完成之后，我们就可以获得被跟踪系统的状
态。cvKa lmancorrect0的结果存储在state_pos仁中，cvKalmanpredict0的结
果存储在sta·e_pre中。
cvKalmanpre clict (
CvKalman · kalman,
const CvMat* control · NULL
cvKalmancorrect(
CvKalman · kalman,
CvMa t"measured
第3页
OCR文字识别结果:
Kalman滤波器示例代码
很明显现在是给出一个例子的时间了。我们举一个相对简单的例子并实现。假设有
一个做圆周运动的点，就像一辆在赛道上行驶的汽车。汽车沿着赛道以几乎恒定的
速度行驶，但是也存在一些波动(也就是过程噪声)。我们使用某种方法测量汽车的
位置，例如用视觉的方法对其进行跟踪。这也会产生一些(不相关的并且可能不同
的)噪声(也就是测量噪声)。
这样，我们的模型就比较简单汽车在任意时刻都具有位置和角速度。这两个元素
合起来形成一个二维的状态向量xk。但是，我们只能测量汽车的位置，所以便得
到一个一维的。向量。zk.
编写一段程序(例10-2)，其输出显示了汽车的圆形运动(红色)，我们的测量值(黄色!
和由Kalman滤波器预测的位置(白色)。
我们在调用的开始包含几个库的头文件。同时，定义一个将被证明是有用的宏，用
于将汽车的位置从角度转换为笛卡儿坐标以便在屏幕上绘图。【359】
倒10-2Ka】..滤波器示例代码
#include'cv. h"
#include'highgui. h"
#include cvx_clefs. h"
#define phi2×y ( mat ) /
cvpoint ( cvRound(img-»width12 lmg-»width13*cos ( mat->
cvRound( img-»heightl2 img-»width13*sin(mat-»data. F1 [0]]
int main{int arge, char** argv )
// Initialize, create Kalman Filter object, window, random number
// generator etc.
cvNamedwindow('Kalman", 1
'"continued below
接下来，我们要创建一个随机数产生器，一幅用于绘制的图像和Kalman滤波器结
构。注意，需要指定Ka】man滤波器状态变量的维数(2)和测量变量的维数(1)。
394 Wiog
第4页
OCR文字识别结果:
.. continued from above
CvRandstate rng
cvRandlnit( & rng, 0, 1, 1, CV_RAND_Ulli )
Iplimage * img-cvcreatelmage« cvsize«500 500 ) 8 3 ) I
/ II
CvKa lman* kalman-cvcreate Kalman( 2, 1, 0 )
. continued below
一旦有了这些块之后，我们为状态x一k、过程噪声讨一k、测量值z_k和最重要的传
递矩阵F创建矩阵(实际上是向量，但在OpencV中我们把所有的都称为矩阵)。状
态需要被初始化，我们将其初始化为一些分布在0附近的合理的随机数。
由于传递矩阵联系着系统在时间七和时间t+1的状态，所以它是至关重要的。在
此例中，传递矩阵是2×2的(因为状态向量是二维的)。实际上，传递矩阵赋予了状
态向量元素意义。我们将x_k视为代表汽车的角度位置((p)和角速度(aJ)。这样，传
递矩阵就为[[1,dt]，[O,川。因此，当乘以F之后，状态(忆co)变成(伊+国dt,a))——
即角速度保持不变而角度位置增加角速度乘以时间步长。在我们的例子中，为了
方便选取dt-1.0，但是实际上需要选取类似连续视频帧间的时间的值。
. cont Inued from above
// state is (phi, delta_phi ) angle and angular velocity
/ / Initialize with random guess.
CvMat* x-k. cvcreate Matl 2, 1, CV-32FC1 )
c'vR andset Rangel & rng, O, 0. 1, 0 )
rng. Di s t type = CV-RAND-NORMAL
cvRand ( & rng, x-k )
/ / process noise
CvMat* w-k · cvcreate Matl 2, 1, CV-32FC1 )
/ / measurements, only one parameter for angle
CvMat* z-k. cvcreate Matl 1, 1, CV-32FC1 )
cvz e ro ( z-k )
// Transition matrix F'describes relationship between
// model parameters at step k and at step k+1 ( this is
/ / the'dynamics"in our model )
跟踪与运动395
第5页
OCR文字识别结果:
const float F[] { 1, 1, 0, 1 }
memcpy ( kalman-»transition_matrix-»data. 61,-F, sizeof ( F ) )
. continued below
[360-36]]
Kalman滤波器还有必须被初始化的其他内部参数。特别地，l×2的测量矩阵H被
一个不直接的恒等函数初始化为【1,0]。过程噪声和测量噪声的协方差被设置为合
理的但是比较有趣的值(你可以自己试一试)，后验误差协方差也被初始化单位矩阵
(这样做是为了保证第一次迭代有意义$它的值随后会被改写)。
类似鹰，由于在此时没有任何信息，我们将(第一个前假象阶段的)后验状态初始化
为一个随机值。
. cont I nued from above
/ / Initialize other Kalman filter parameters.
cvsetldentity( kalman-»measurement_matrix,
cvRea1Scalar ( 1 ) )
cvsetldentity( kalman-»process_noise_cov, cvRea1Scalar ( le-5 ) )
cvsetldentity( kalman-»measurement_noise_cov, cvR. ea1scalar ( le-
1 ) )
cvsetldentity( kalman->error_cov_post, cvRealscalar ( 1 ) )
while( )
. continued below
最后，我们可以在实际动态系统上开始预测了。首先，用Kalman滤波器预测它所
认为在此阶段会产生的什么(也就是在获得任何新的信息之前);我(门称为y_k。然
后为这次迭代产生z_k(测量值)的新值。根据定义，这个值是。实际的。值x_k乘
以测量矩阵H再加上随机测量噪声。必须在这里提一下，除了类似这个简单应用
之外，不能从x一k产生z_k，而是由环境状态或传感器来产生一个生成函数。在这
个模拟的例子中，我f门添加随机噪声从一个基本的。实际。数据模型中产生测量
值这样就可以看到Kalman滤波器的作用了。【361～3621
. continued from above
/ / predict point position
const CvMat"y-k-cvKalmanpredict( kalman, O )
396 8108
'勇
第6页
OCR文字识别结果:
/ / generate measurement ( z-k )
cvRandset Range (
& rng,
sqrt ( kalman->measurement-noise-cov->data. F1 [0] ),
cvRand( & rng, z_k )
cvMatMu1Add( kalman-»measurement_matrix, x_k, z_k, z_k )
. continued below
画出由先前合成的观测值、由Kalman滤波器预泖」的位置和基本状态(在这个模拟的
例子中我们恰巧知道)所代表的三个点。
. cont inued f rom above
// plot points ( eg convert to plainar coordinates and draw)
cvzero( img )
cvcircle( img, phi2 ×y(z-k ), 4,
/ / observed state
cvcircle ( · img, phi2 ×y (y-k ), 4,
/ /"predicted"state
cvcircle ( img, phi2×y ( x-k ), 4,
cvshowlmage ("Kalman", img )
. continued below
CVX-YELLOW )
CVX_WHITE, 2 )
CVX-RED ) I / / real state
到此，我们可以开始下一次迭代了。首先要做的是在此调用Kalman滤波器并赋予
其最新的测量值。接下来就是产生过程噪声。然后对x_k乘以时间传递矩阵F完
成一次迭代并加上我们产生的过程噪声现在，我们又可以开始新的一轮计算。
. continued from above
/ / adjust Kalman filter state
cvKalmancorrect( kalman, z-k )
// Apply the transition matrix'F'( e. g., step time forward )
/ / and also apply the"process"noise w k
cvR. ands e t R. ang e
& rng,
跟踪与运动397
第7页
OCR文字识别结果:
sqrt ( kalman->process-noise-cov->data. F1 [0] ),
◆
cvRand( & rng, w_k )
cvMatMu1Add ( kalman-»transition_matrix, x_k, w_k, x_k )
/ / exit if user hits Esc'
if ( cvwait Key( 100 ) 27 ) break
return o
[362-363]
可以看出，Kalman滤波部分并不十分复杂，有一半所需的代码只是用来产生一些
填充信息。在任何情况下，我们都应该总结做过的所有事情，为的是确保所有的都
有意义。
首先创建表示系统状态的矩阵和将要获得的测量值的矩阵开始，定义传递矩阵和测
量矩阵，然后初始化噪声协方差矩阵和滤波器的其他参数。
在初始化状态向量为一随机值后，调用Kalman滤波器并让其做出第一次预测。一
旦获得此预测值(虽然此值在第一次中没有太大的意义)，将其在屏幕上绘制出来。
同时也合成一个新的观测并在屏幕上绘制出来与滤波器的预测值进行比较。然后，
以新测量值的形式传递新信息给滤波器，而这个新的测量值会被合并到滤波器的内
部模型中。最后，给模型合成一个新的。真实的。状态，这样我们就可以再次重复
这个循环。
运行代码，可以看到小的红色球绕轨道一圈一圈地运动。小的黄色球围绕红色球出
现和消失，这表示Kalman滤波器正尝试去。看透。的噪声。白色的球收敛在红色
球周围的小空间中运动，这说明Kalman滤波器在我们的模型中对物体(汽车)的运
动做出了合理的估计。
在例子中还没有讨论的一个问题是控制输入的使用。例如，如果这辆车是一辆无线
控制的汽车，同时我们知道人会控制其进行什么操作，那么就可以将这个信息加入
模型中。这样，速度就由控制器来设定。我们需要补充一个矩阵B(kalman一
〉con七ro1_mat八·)，同时还要为cvKalmanpre己ict0调节控制向量u提供第二个
参数。
398第lo章
一型勇
第8页
OCR文字识别结果:
肾匡
扩展Kalman滤波器简述
你可能会注意到，要求系统的动力学性能在参数上呈线性是相当苛刻的。当动力学
性能不是线性时，Kalman滤波器对我们仍然有用，OpencV中的Kalman滤波器函
再说一T,。线性。的意思(实际上)是Kalman滤波器定义中的不同的阶段可以用
矩阵来表示。什么时候情况又不是这样的呢?实际有许多可能。例如，假设我们的
控制测量是汽车的油门被踩下次数的总数汽车速度和油门被踩的次数之间的关系
就不是线性关系。另一个比较常见的问题是施加在汽车上的力容易用笛卡儿坐标系
来表示，而汽车的运动(如我们举的这个例子)更容易用极坐标来表示。这个问题可
能在这样的情况中出现，即当用船代替汽车在规则的水流中做圆形运动但同时又向
着某一特定的方向前进时。【363-364】
在所有这样的情况中，仅仅使用Kalman滤波器是不够的。处理这些非线性(或者试
图去处理)的一种方法是将相关的过程(例如，更新F或者控制输人响应B)线性化。
这样，我们需要根据状态x在每一个时1司步长计算F和B的新值。这些值只是在
某个特定x值的附近接近真实的更新和控制函数，但是在实际中已经足够。这个
对Kalman滤波器的扩展被简单地叫作扩展Kalman滤波器[Schmidt66]。
OpencV没有提供专门的程序来实现扩展Kalman滤波器，但实际上我们不需要这
样一个函数。我们所需要做的只是在每次更新前重新计算和设置kalman一
》up己a上e_matrix和kalman-〉con仁ro1_matrix的值。Kalman滤波器后来由一个
unscented粒子滤波[Merwe00]的公式被巧妙地扩展到非线性系统。[Thrun05]非常清
楚地介绍了Kalman滤波器的全部信息，包括最新的发展。
condensation 算法
Kalman滤波器是基于单个假设的建模。由于此假设的概率分布的基本模型是单高
斯的，所以不可能用Kalman滤波器同时表达多个假设。要解决这一问题，可以使
用一个更先进一点的方法，即condensation算法。condensation算法是建立在称为
粒子滤波器的一类更广泛的预估器的基础上的。
为了理解condensation算法的意图，假想有一个物体以恒定的速度运动(如Kalmari
墓设。现在假设有一个被遮挡的运动的物体。我们并不清楚这个物体怎样在运
跟踪与运动399
良医
第9页
OCR文字识别结果:
动，它可能以恒定的速度继续运动，可能停止然后/或者改变运动方向。简单地增
大物体位置(高斯)分布的不确定性并不能使Kalman滤波器描述这些多种可能性。
Kalman滤一波器必须为高斯的，所以它不能表示这种多态分布。
和Kalman滤波器相同，有两个函数(分别)用于创建和销毁表示condensation滤波
器的数据结构。惟一的不同在于对于condensation滤波器，创建函数cvcreate一
con Densationo多了一个参数。输给参数的值为滤波器在任意给定时刻都保持不
变的假设(也就是。粒子。)的数量。这个数字应该相对比较大(50或者100，对于
复杂的情况可能更大)，因为这些单个的假设将替代Kalman滤波器的参数化高斯概
率分布。如图10-20所示。【364～365】
Cvcon Densat ion* cvcreatecon Densat ion (
int dynam-params,
int measure-params,
int sample-count
void cvReleasecon Densation(
Cvcon Densation"condens
(a)(匕)
图1 0-20可以(a)和不可以(b)用均值和协方差参数化的连续高斯分布表示的分
布两个分布都可以一组粒子来代替，这组粒子用密度来近似概率分布
这个数据结构含有如下内部元素
typedef struct Cvcon Densation
int MP // Dimension of measurement vector
int DP // Dimension of state vector
float"DynamMatr // Matrix of the linear Dynamics system
float"State // Vector of State
400第10章
月
第10页
OCR文字识别结果:
P
in 匕 SamplesNum
float** f1Samples
£1oat** £1Newsamples
£1oat* f1Conf idence
float* £1Cumulativel
float* Temp
f 1 oa t · Randoms amp 1 e
CvRandstate* Rands
/ / Number of Samples
/ / array of the Sample Vectors,
/ / temporary array of the Sample
/ / Vectors
/ / Con£idence for each Sample
/ / Cumulative con£idence
/ / Temporary vector
/ / Randomvector to update sample set
/ / Array of structures to generate
/ / random vectors
Cvcon Dens ation
分配了condensation滤波器数据结构后需要初始化这个结构。我们用
cvcon Denslnisampleseto来实现。在创建cvcon Densation结构时，需要指出
粒子的数目并为每个粒子指定维数。初始化所有的粒子是相当麻烦的®。幸运的
是，cvConDenslnitsampleset 0以一种方便的方法完成了这个任务我们只需要
为每一个维度指定一个范围。
void cvcon Denslnitsampleset(
Cv Con Densation+ condens,
CvMat · lower-bound,
CvMat · upper-bound
这个函数需要我们初始化两个cvMat结构。这两个结构均是向量(意思是它们都只
有一列)，并且每个都有与系统状态维数相同的元素个数。这两个向量将被用于设
置初始化cvcon Densation结构中粒子向量的范围。
下面的代码创建了两个大小为Dim的矩阵并将它们分别初始化为-1和+】。当调用
cvcon Densln上匕sampleset o时，初始的粒子每个将被初始化为在-1到+】1司0的
随机数。所以，如果Dim为3，则滤波器被初始化为均匀分布在中心在原点、边长
为2的立方体内部的粒子。
CvMat LB · cvMat{Dim, 1, CV-MAT32F, NULL )
CvMat UB. cvMat(Dim, 1, CV-MAT32F, NULL )
cvmA11oc ( &LB )
cvmA11oc ( &UB )
当然，如果你了解粒子滤波，那么你就知道这里正是我们用有关系统状态的先验知识
(或先验假设)初始化滤波器的地方。初始化滤波器的这个函数只是帮助我们产生一个
点集的均匀分布(也就是概率密度函数是平的)。
跟踪与运动401
第11页
OCR文字识别结果:
cvcon Denslnitsampleset ( Con Dens, & LB, &UB )
最后，一个用于更新condensation滤波器状态的函数
void cvcon DensuipdateByTime( Cvcon Densation · condens )
使用这个函数还需要比我们表面所看到的多做些工作。具体说来，必须根据前一次
更新之后可获得的信息来更新所有粒子的置信度。不好的是，OpencV中没有为此
而设置的方便的函数。原因是粒子新的置信度和新的信息间的关系取决于应用的环
境。这里有一个更新的例子，它对滤波器中的每个粒子的置信度进行了简单更
// Update the con£idences on all of the particles in the filter
// based on a new measurement Mt]. Here M has the
// dimensionality of
// the particles in the filter.
void CondprobDens(
Cvcon Densation* CD,
float* M
£or( int I-0 I<CD-»SamplesNum i++
float p 1. 0f
£or( int j = 0 j«CD->DP j++ )
p ( £1oat ) expl
-0. 05 · (M [ j ] ー CD->£ 1Samples [ I ] [ j ] ) · (M [ j ] ー CD->
f1samples[i] [i]]
P*i
CD-> f ICon f idene e [ 5L ] Prob
®细心的读者会注意到过个更新实际上暗含了一个高斯概率分布，但是你当然可以根据
你专门的应用环境选用一个更为复杂的更新。
402第10章
第12页
OCR文字识别结果:
当置信度被更新后就可以调用cvcondensupdate By Time O去更新粒子了。这里
。更新。的意思是重采样，即产生与计算出的置信度一致的一组新的粒子。更新之
后，所有的置信度又将变成1.0f，但是粒子分布将会把前一次修改的置信度直接包
含进下一次迭代粒子的密度中。【366～367】
练习
在.../opencvTsamptes/c/目录中是一些示例代码程序，演示了本章讨论的很多算法
lkdemo. c l%m
cansa垆demo.c(用mean-shift根据颜色跟踪某个区域)
motempl.c(运动模板)
*kalman.c(Kalman滤波器)
1.cvGood FeaturesTo Track0中的协方差Hessian矩阵是在图像中的正方形区域
上计算得出的正方形区域大小由此函数的block_si·e设置。
a、从概念上来讲，当block_size增加是会出现什么情况?我们会得到更多
的还是更少的。好的特征。?为什么?
b，仔细阅读lkdemo.c代码，找到cvGood FeaturesToTracko，改变
block_size的大小，然后运行查看有何不同。
2.参考图10-12，思考实现寻找亚像素级角点的函数cvFindcornersubpixo。
a，在图10-12中，如果棋盘格扭曲了，直的黑白线形成交于一点的曲线时，
会发生什么情况?亚像素级角点还能找到吗?请给出解释。
b，如果扩大扭曲的棋盘格的角点周围窗口的尺寸(扩大win和zero_zone参
数)，寻找亚像素会更精确呢还是它会开始发散?请解释你的答案。
3光流
a，描述一个用块匹配跟踪比Lucas-Kanade光流跟踪更好的物体。
b，描述一个用Lucas-Kanade光流跟踪比块匹配更好的物体。
4.编译tkdemo.c.连上一个网络摄像机(或者用以前拍摄的一个有纹理的运动物
体的序列)。在运行此程序中，。，。是自动初始化跟踪，。c。是清除跟踪，
单击鼠标添加一个新的点或删除一个已有的点。运行lkdemo.c,同时键人
跟踪与运动403
'良
第13页
OCR文字识别结果:
。r。初始化点的跟踪。观察效果。
a，现在进入代码，去掉亚像素级角点定位函数cvFindcornersubpix()。这
样做影响结果吗?怎样影响的?
b，再次进入代码，不调用cvGood FeaturesTo Tracko，在物体周围的ROI
中标记一些格子的点。描述这些点会产生什么变化并说明为什么。
提示所发生的部分是由于孔径问题造成的——给定一个固定窗口尺寸和
一条直线，我们分辨不出这条线是如何运动的。
5 6 7
修改tkdew工程序，使其对轻微的摄像机运动进行的简单图像稳定。在比摄
像机的输出大得多的窗口中心显示稳定的结果(当第一组点保持稳定时帧画面
可能抖动)。
用一个网络摄像机或一个运动的有颜色的物体的彩色视频，编译和运行
camsh垆de科o.c.用鼠标画一个(刚好)包含这个运动物体的方框这个程序将会
跟踪此物体。
a，在eamsh垆dento.c中，用cvMeanshift0替换cvcamshifto函数。描述
一个跟踪器比另一个跟踪器跟踪效果更好的情况。
b，编写一个在初始的cvMeanshift 0区域中标记一组点的函数。立即运行两
个跟踪器。
c.这两个跟踪器如何结合使用会让跟踪更稳定?解释原因并且演示。
使用网络摄像机或者事先存储的视频文件，编译并运行运动模板代码
motempl. c.
a，修改motempi.c,让它可以进行简单的姿态识别。
b，如果摄像机在运动，阐述一下怎样使用练习5中的运动稳定的代码让运动
模板对轻微运动的摄像机也有效。
8.描述一下用一个线性状态模型(非扩展)Kalman滤波器怎样去跟踪圆周(非线性)
运动。
提示怎样预处理这种情况来得到线性的动态?
9.
建立一个当前状态取决于先前状态的位置和速度的运动模型。将tkde科o.c(只
使用一些click点)结合Kalman滤波器更好地跟踪Lucas-Kanade点。显示每一一
个点的不确定性。这个跟踪在哪里会失败?
提示将Lucas-Kanade作为Kalman滤波器的观测模型，调整噪声使其可以跟
404第10章
。..。-：.一。一
第14页
OCR文字识别结果:
P
-
踪。保持运动的合理性。
10.
Kalman滤波器依赖于线性动态性和Markov独立性(也就是它假设当前状态只
依赖于刚过去的状态而不是所有过去的状态)。假设要跟踪一个运动与先前的
位置和速度相关的物体，但是你错误地使用了一个动态项使状态依赖于先前的
位置——换句话说，遗漏了先前的速度项。
a.Kalman的假设仍然还成立吗?如果成立，说明为什么?如果不成立，说明
为什么不满足假设。
b，怎样建立Kalman滤波器使其在动态项的某些部分被漏掉时仍可以进行
跟踪?
提示想一想噪声模型。
11.使用一个网络摄像机或者一段视频，内容为一个人挥舞双手，每个手中都有一
个明亮颜色的物体。用condensation去跟踪两只手。
跟踪与运动405