第六章
光流分析
6.1 运动场 (Velocity)与光流 (Optical Flow)
一 .运动场
物体与摄象机之间的任何相对运动都将导致视平面上与空间物体对应点发生变化。严格地说空
间运动在视平面上引起的对应运动称为运动场,如图 6 .1 所示。假定空间点 P 与摄象机之间相
对运动的速度为 V,在视平面上对应点 p 的运动速度为 v,则
dt
drv
dt
dRV ==
(6.1)
且由投影关系有
),( ZR
Rfr r=
(6.2)
从而在投影关系下有
22 ),(
)(
),(
),(),(
ZR
ZVR
ZR
RZVVZR
f
v r rr rr ××=?=
(6.3)
O o
X,Y x,y
Z
P
p
R
r
v tδ
V tδ
图 6 .1 空间运动及其运动场
二 .光流
在图象中可测的是图象辐照度的变化,这种变化并不总是与物体和摄象机间的
相对运动对应的。可以从下面的例子中看出这一点。典型的具有均匀反射特性的球
绕其一轴旋转,从摄得的图象上看不出任何变化;反过来光照条件的变化却会引起
图象的变化,如图 6 .2 所示。
因此这种从图象上测得的变化并不完全反映实际的运动情况,称之为表观运动
(Apparent motion)。所谓光流是指亮度模式引起的表观运动,理想的情况是这种表观
运动反映了实际的运动。
图 6 .2 表观运动
6.2 基于梯度的光流计算
在基于梯度的方法中,时空梯度之间的关系是极其重要的,这个关系被称之为基本
等式,它构成了对光流计算的一个重要约束。
6.2.1 基本等式的引入
设图象上的点 ( , )x y 在时刻 t 的辐照度为 f x y t( , , ) ,经过间隔 ?t 后对应点为
f x x y y t t( , , )+ + +? ? ? ,当 ?t → 0时可以认为辐照度不变,于是有
f x y t f x x y y t t( , , ) ( , , )= + + +? ? ? (6.4)
由 Taylor 展开
f x x y y t t f x y t fx x fy y ft t( , , ) ( , , )+ + + = + + + +? ? ? ? ? ??? ?? ?? ε
(6.5)
忽略二阶无穷小,由于 ?t → 0,于是
?
?
?
?
?
?
f
x
dx
dt
f
y
dy
dt
f
t+ + = 0
f u f v fx y t+ + = 0
(6.6)
其中 u dxdt v dydt= = ,分别是在 x 和 y 方向上的光流的分量
对基本等式的理解
理解之一 :
式 (6.6)可理解为
df
dt = 0。
理解之二 :( )
( ) tyx fvuff ?=? ,于是沿 ( )yx ff 方向
上光流的大小为
f
f f
t
x y
2 2+ ;
理解之三 :
将 (6.6)式看成是 uv 平面上的一条直线,一个方程
两个未知数,无法定解。为了定解必须附加上其
他约束。
v
u
f , fx y
约 束 线
图 6 .3 基本等式所确定的约束线
6 .2.2 光流计算与附加约束
基本等式给出了光流计算的一个约束,但仅有这一个方程,是无法确定两个未
知量 u 和 v 的,必须引入其它附加的约束才有可能唯一确定光流场。研究者们从不
同的角度出发引入不同的约束从而导致产生不同的光流分析方法,下面讨论几种典
型的方法。
6.2.2.1 基于一阶梯度的方法
Horn 和 Schunck 所采用方法基本思想是在求解光流时,光流本身尽可能平滑,
即引入对光流的平滑性约束。设平滑性约束项为
( )∫∫ +++= dxdyvvuuE yxyxs 2222 (6.7)
由基本等式 (6.6),显然要求
( )∫∫ ++= dxdyfvfufE tyxc 2 (6.8)
于是由 (6.7)和 (6.8)可知,最后求得的光流应满足 (6.9)式,即
( )[ ]{ }∫∫ ++++++= dxdyfvfufvvuuE tyxyxyxs 22222min l (6.9)
对形如 (6.10)形式的形式变分问题{ }
dydxvvuuvuF yxyx∫∫ ),,,,,(min (6.10)
的解是对应 Euler 方程 (6.11)的解
?
?
?
???
?
=??
=??
0
0
y
F
x
FF
y
F
x
FF
yx
yx
vv
v
uu
u
?
?
?
?
?
?
?
?
(6.11)
对于 (6.9)式
( )22222 tyxyxyx fvfufvvuuF ++++++= l (6.12)
于是对应的 Euler 方程为
??
???
++=?
++=?
)(
)(
2
2
tyxy
tyxx
fvfuffv
fvfuffu
l
l (6.13)
其中 ?
2
是 Laplace 算子。
(6.9)式中 λ反映对数据及约束的信度。当数据本身含有较多的噪声时,则原始
数据的可信度较低,更多地依赖对光滑性的约束, λ可以取较小的值,反之 λ可以取
较大的值。
对于 Euler 方程 (6.13)进行研究不难发现对于以下几种情况难以定解:
l 当区域亮度梯度为零时,无法确定光流;
l 当物体沿某一边缘运动时,边缘上的光流无法确定;
l 图象的四周、角点梯度变化较快,计算的光流值有较大的偏差。
解决的方法 -- 内插
实际计算时由于计算的对象是离散化的图象,因此需要进行离散化处理。
离散化后光滑性约束变成
s u u u u v v v vij i j i j i j i j i j i j i j i j= ? + ? + ? + ?+ + + +14 1 2 1 2 1 2 1 2( ) ( ) ( ) ( ), , , , , , , ,
而基本等式的约束变成 ( )
2
tijyijxij fvfufc ++=
于是极小化目标函数为
?
?
?
?
?
? += ∑∑
i j
ijij cse )(min l (6.14)
对 (6.14)求关于 ukl和 vkl的偏导,并令其为零,整理后有
??
???
?=++
?=++
tyklklyklxy
txklklyxklx
ffvvfuff
ffuvffuf
lll
lll
)1(
)1(
2
2
(6.15)
其中 ukl、 vkl分别是 ukl和 vkl的四邻域平均
解 (6.15)得
?
?
?
???
?
++
++?=
++
++?=
y
yx
tklyklx
klkl
x
yx
tklyklx
klkl
fff fvfufvv
fff fvfufuu
)(1
)(1
22
22
l
l
于是得到自然的迭代过程
?
?
?
???
?
++
++?=
++
++?=
+
+
y
yx
t
n
kly
n
klxn
kl
n
kl
x
yx
t
n
kly
n
klxn
kl
n
kl
fff fvfufvv
fff fvfufuu
)(1
)(1
22
1
22
1
l
l (6.16)
上述迭代过程如图 6 .4 所示。
v
u
f , fx y
约 束 线 u , v( )n n
u , v +1n+1n( )
图 6 .4 迭代求解过程
6.2.2.2 基于高阶梯度的方法
与 Horn 等的方法不同, Tretiak 等考虑图象本身在灰度上连续性,对灰度场加约束。
由基本等式有
?
?
?
?
?
?
f x dx y dy t
x u x dx y dy t
f x dx y dy t
y v x dx y dy t
f x dx y dy t
t
( , , ) ( , , ) ( , , ) ( , , )
( , , )
+ + + + + + + + + +
+ + + = 0
(6.17)
对 (6.17)式中各项在 ( , , )x y t 点利用 Taloy 展开,有
?
?
?
?
?
?
?
? ?
f x dx y dy t
x
f x y t
x
f x y t
x dx
f x y t
x y dy
( , , ) ( , , ) ( , , ) ( , , )+ + = + +2
2
2
?
?
?
?
?
? ?
?
?
f x dx y dy t
y
f x y t
y
f x y t
y x dx
f x y t
y dy
( , , ) ( , , ) ( , , ) ( , , )+ + = + +2 2
2
?
?
?
?
?
? ?
?
? ?
f x dx y dy t
t
f x y t
t
f x y t
t x dx
f x y t
t y dy
( , , ) ( , , ) ( , , ) ( , , )+ + = + +2 2
u x dx y dy t u x y t u dx u dy
v x dx y dy t v x y t v dx v dy
x y
x y
( , , ) ( , , )
( , , ) ( , , )
+ + = + +
+ + = + +
于是 (6.17)等价于
( ) ( )
( ) ( )
( ) ( ) 02
2
=++++++
++++++++
++++++++
dyvfufdxdyvfvfufuf
dxvfufdyfvfufvfuf
dxfvfufvfuffvfuf
yyyyxyyxyxyyyxxxxy
xyxxxxtyyyyxyyxy
txxyxxyxxxtyx
(6.18)
由基本等式及 dx dy, 取值的任意性于是有
f u f v f
f u f v f u f v f
f u f v f u f v f
f u f v
f u f u f v f v
f u f v
x y t
xx yx x x y x tx
xy yy x y y y ty
xx x yx x
xy x xx y yy x xy y
xy y yy y
+ + =
+ + + + =
+ + + + =
+ =
+ + + =
+ =
0
0
0
0
0
0 (6.19)
上述六个方程并不独立,直接求解有困难。由于光流场本身的在大多数区域是连续
平滑的,因此可以假定 u u v vx y x y, , , 近似为零,则有
f u f v f
f u f v f
f u f v f
x y t
xx yx tx
xy yy ty
+ + =
+ + =
+ + =
0
0
0 (6.20)
在 (6.20)式中两个未知数三个方程,可以用最小二乘法进行计算。
6.3 基于区域匹配的光流计算方法
与基于梯度的方法不同,在这种方法中考虑了局部保守 (Conservation)信息,即局部上的点
具有相似的速度。利用这一信息得到对光流的初始估计,并可以根据这一估计值利用邻域
关系得到最终的速度场。
一 . 基于保守信息的初始速度场估计
设相邻的两帧图象分别为 f i j1( , ) 和 f i j2( , ) ,考虑第一帧图象上的点 f x y1( , ) 周围大小为
( ) ( )2 1 2 1n n+ × + 邻域,考虑对应的运动向量为 ( , )u v ,且 ? ≤ ≤N u v N, 。实际上 N 确定了
搜索的范围,可以用 (6.23)式衡量对应的误差 E u vc ( , ) 。
E u v f x i y j f x u i y v jc
i j n
n
( , ) ( , ) ( , )
,
= + + ? + + + +
= ?
∑ 1 2
2
(6.23)
小 小
小大 大
大
(a) 变化平缓部分 (b) 边缘附近 (c) 角点附近
图 6 .11 不同类型区域对应的运动向量的信度
上式的误差可以转化转化为对所有可能运动向量 ( , )u v 的分布的描述可以定义对于运动向量为
( , )u v 的权为 ,
W u v ec kE u vc( , ) ( , )= ? (6.24)
于是可以得到加权的运动向量
u
W u v u
W u v
v
W u v v
W u v
c
c
u N
N
u N
N
c
v N
N
u N
N c
c
u N
N
u N
N
c
v N
N
u N
N= =
= ?= ?
=?= ?
= ?= ?
=?= ?
∑∑
∑∑
∑∑
∑∑
( , )
( , )
,
( , )
( , )
(6.25)
事实上,由于某些区域得到的运动向量并非是可靠的,例如对于平坦的区域就很难得到准确的
运动向量。图 6 .11 是对于不同类型区域的运动向量分布情况,其中正交的直线是对应速度场
协方差矩阵 (4.26)式的特征向量方向。
( ) ( )( )
( )( ) ( )
??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
???
???
=
∑ ∑
∑ ∑
∑ ∑
∑ ∑
∑ ∑
∑ ∑
∑ ∑
∑ ∑
?= ?=
?= ?=
?= ?=
?= ?=
?= ?=
?= ?=
?= ?=
?= ?=
N
Nu
N
Nv
c
N
Nu
N
Nv
cc
N
Nu
N
Nv
c
N
Nu
N
Nv
ccc
N
Nu
N
Nv
c
N
Nu
N
Nv
ccc
N
Nu
N
Nv
c
N
Nu
N
Nv
cc
c
vuW
vvvuW
vuW
vvuuvuW
vuW
vvuuvuW
vuW
uuvuW
C
),(
),(
),(
),(
),(
),(
),(
),(
2
2
(6.26)
对于变化平缓的区域,与特征向量对应的特征值均很小,而在角点附近对应的特征值均较
大,实际上特征值的大小反映了对速度场估计值的信度。因此有必要利用邻域信息对运动向量
进行求精。
二 . 利用邻域信息的求精
运动点周围邻点的速度场对当前点速度值提供了一种支持,但这种支持不是等权的。实际
上可以考虑两种决定因素,即周围邻点的距离和邻点本身的信度,前者可以利用 Gauss 函数作
为加权,后者可以考虑利用从 (6.26)式的协方差阵中得到的特征值进行加权。
取大小为 ( ) ( )2 1 2 1w w+ × + 邻域,利用上述原则可以得到权函数 W i jn ( , ) ( , )? ≤ ≤w i j w ,
于是有利用邻域信息的当前点加权速度值
u
W i j u
W i j
v
W i j v
W i j
n i j
j w
w
i w
w
n
j w
w
i w
w
n i j
j w
w
i w
w
n
j w
w
i w
w= =
=?=?
=?=?
=?=?
=?=?
∑∑
∑∑
∑∑
∑∑
( , )
( , )
,
( , )
( , )
, ,
(6.27)
其中 ( , ), ,u vi j i j 为对应点 ( , )i j 的理想速度值,并设当前点的坐标为 (0,0)。
对应的协方差阵为
( ) ( )( )
( )( ) ( )
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
=
∑ ∑
∑ ∑
∑ ∑
∑ ∑
∑ ∑
∑ ∑
∑ ∑
∑ ∑
?= ?=
?= ?=
?= ?=
?= ?=
?= ?=
?= ?=
?= ?=
?= ?=
w
wi
w
wj
n
w
wi
w
wj
jin
w
wi
w
wj
n
w
wi
w
wj
jijin
w
wi
w
wj
n
w
wi
w
wj
jijin
w
wi
w
wj
n
w
wi
w
wj
jin
n
jiW
vvjiW
jiW
vvuujiW
jiW
vvuujiW
jiW
uujiW
C
),(
),(
),(
),(
),(
),(
),(
),(
2
,,,
,,
2
,
(6.28)
现在从式 (6.25)和 (6.27)可以得到两个速度场的估计值,考虑最优估计 U u v= ( , ) 应满足
与上述两个估计值的距离尽可能小,即
( ) ( ) ( ) ( )[ ]{ }∫∫ ??+?? ?? dxdyUUCUUUUCUU ccTcnT 11min (6.29)
利用变分有
C U U C U Un c c? ?? + ? =1 1 0( ) ( ) (6.30)
事实上由于无法事先得到理想的速度值,因而实际上不能通过 (6.27)和 (6.28)式计算出
U 和 Cn ,只能通过松弛方法进行计算。利用 Gauss-Siedel 迭代有
U C C C U C U
U U
k c n c c n k
c
+
? ? ? ? ?= + +
=
1
1 1 1 1 1
0
(6.31)
6.4 从运动恢复结构
6.4.1 空间运动与对应视平面的运动场
在讨论空间运动时我们假定摄象机和外界环境 (这里的外界环境是指摄象机中的景物 )
只能有一个是运动的,不失一般性我们假定摄象机是运动的,与第三章类似假定摄象机
仅按 ( )TCBAh = 平移,按 R 旋转。设 t 时刻空间点 F X Y Z t( , , , ) 经过时间 ?t 后在摄
象机坐标系下为 F X X Y Y Z Z t t( , , , )+ + + +? ? ? ? 。则由 (3.49)式有
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=
?
?
?
?
?
?
?
?
?
?
?+
?+
?+
C
B
A
Z
Y
X
R
ZZ
YY
XX
T (6.46)
由 (3.44)式知
?
?
?
?
?
?
?
?
?
?
??+??+??????
??????+??+??
?+????????+?
=
)cos1(cossin)cos1(sin)cos1(
sin)cos1()cos1(cossin)cos1(
sin)cos1(sin)cos1()cos1(cos
2
3123213
132
2
2312
231321
2
1
lllllll
lllllll
lllllll
R
在 0 点处对 ?作一阶 Taloy 展开,有 :
)(
1
1
1
2
12
13
23
?+
?
?
?
?
?
?
?
?
?
?
?????
?????
?????
= O
ll
ll
ll
R (6.47)
由于 ?? ?= =ω ω ωt li i, ,于是有
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+=
?
?
?
?
?
?
?
?
?
?
?+
?+
?+
C
B
A
Z
Y
X
tI
ZZ
YY
XX
0
0
0
12
13
23
ww
ww
ww
考虑投影关系,取 ?t → 0 有
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
×
?
?
?
?
?
?
?
?
?
?
?=
?
?
?
?
?
?
?
?
?
?
c
b
a
f
y
x
f
Z
Z
Y
X
3
2
1
w
w
w
&
&
&
(6.48)
为书写方便,将
?
?
?
?
?
?
?
?
?
?
dtdZ
dtdY
dtdX 简记为
?
?
?
?
?
?
?
?
?
?
Z
Y
X
&
&
& ,
?
?
?
?
?
?
?
?
?
?
3
2
1
w
w
w 为摄象机旋转的角速度,
?
?
?
?
?
?
?
?
?
?
=
?
?
?
?
?
?
?
?
?
?
dtdC
dtdB
dtdA
c
b
a 为平移速度。
又由于
由投影关系有
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
?
=
?
?
?
?
?
?
?
?
?
?
=
?
?
?
?
?
?
?
?
?
?
0
y
x
f
Z
f
y
x
f
Z
Z
Y
X
dt
d
Z
Y
X
&
&&
&
&
&
(6.49)
为与前面记号上统一,记 u x v y= =&, & ,由 (6.48)和 (6.49)有
???
?
???
? +????
?
??
?
? ?+?=
???
?
???
? +????
?
??
?
? +??=
f
x
f
y
Z
cyxf
Z
bfv
f
x
f
y
Z
cxyf
Z
afu
21
31
21
32
wwww
wwww
(6.50)
4.4.2 平面对象的光流
一 . 平面光流参数
如果在视平面上成象的对象是平面,考虑其方程为 Z pX qY r= + + ,则式 (6.50)变成
( )
( )yFyExfDyCxfVv
xFyExfByAxfUu
++++=
++++=
1
1
(6.51)
其中
U ar V br A pa cZ B qar
C pbZ D qb cr E pcr F qcr
= ? ? = ? + = + = +
= ? = + = ? ? = ? +
ω ω ω
ω ω ω
2 1 3
3 2 1
称为平面光流参数。
在式 (4.51)中,有 8 个变量,当有 4 对以上视平面上的点的位置和光流值已知时可以唯一地
确定平面光流参数和摄象机的运动。
二 . N矢量与速度变换矩阵
由于 N 矢量表示具有其明显的优点,下面考虑以 N 矢量形式表示的运动关系。
由第三章有 ′ =m rT mT ,当 T I= 时, ′ =m m ,对于瞬时微小运动 T 可以写成有
T I W t= + ? ,下面首先确定 &m dm dt= 与 W 的关系。
由于 T 保留有一个常数的比例因子,故不妨设 det[ ]T = 1,从而
det[ ] Tr [ ] (( ) )I W t W t O t+ = + +? ? ?1 2 (6.52)
于是 Tr[ ]W = 0,这样有
′ = +m r I W t m( )? (6.53)
′ = + + =m r r m W m t O tT2 2 2 22 1( , ) (( ) )? ? (6.54)
由 (6.54)式并考虑按等比级数展开,有
r O tm W m t O t m W m t O t
m W m t O t
T
T
T
2
2
2 2
2
1
1 2 1 1 2
1 2
= ++ = + ? +
= ? +
(( ) )
( , ) (( ) ) ( , ) (( ) )
( , ) (( ) )
?
? ? ? ?
? ? (6.55)
由于 (6.55)式中 O t(( ) )? 2 可以忽略因此
r m W m tT= ?1 ( , )? (6.56)
从而 (6.53)式可以写成
[ ]( )
[ ] ))((),(
),(1
2tOtmmWmmWm
mtWItmWmm
TT
TT
?+??+=
?+??=′
(6.57)
于是视平面上的光流采用 N 矢量形式可以表达为
& lim ( , )m m mt W m m W m m
t
T T= ′ ? = ?
→? ?0 (6.58)
并且满足 Tr[ ]W = 0。
W 称为速度变换矩阵。
三 . 平面光流参数与速度变换矩阵的关系
下面讨论平面光流参数 { }FEDCBAVU ,,,,,,, 与速度变换矩阵 W 的关系。
由于 T I W t= + ? ,由第三章有
( )
))(()(1)(
)(1
)(
)1(
)1(
2
231321331131
332313
312111
332313
312111
tOtxyWxWfyWxWWfWx
tWfyWxW
tfWyWxWx
ftWtyWtxW
tfWtyWxtWfx
?+??
?
?
??
? +?+?++=
?+++
?+++=
?++?+?
?+?+?+=′
即
u fW W W x W y f W x W y x= + ? + ? +31 11 33 21 13 231( ) ( )
(6.59)
将 (6.59)式与 (6.51)式对比,有
U W A W W B W E W F W= = ? = = ? = ?31 11 33 21 13 23 (6.60)
类似地处理 v 分量,有
V W C W D W W E W F W= = = ? = ? = ?32 12 22 33 13 23 (6.61)
由 (6.60)和 (6.61)并考虑 Tr[ ]W = 0 可以导出平面光流参数 { }FEDCBAVU ,,,,,,, 与速
度变换矩阵 W 之间满足
( )
?
?
?
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
??+
?=
?
?
?
?
?
?
?
?
?
?
+?
?+?
??
=
0
0
0
1
1
100
010
001
3
3)(
3)2(
3)2(
12
13
23
ww
ww
ww
cbaq
p
rr
cqbpa
DAVU
FDAB
ECDA
W
(6.62)
四 . 速度变换矩阵 W的计算
讨论利用 N矢量直接计算速度变换矩阵 W,从 (6.58)式知速度变换矩阵 W应满足
??
???
??
??? ?+?∫
? )(
2 )(),(min
S
TT mdmmWmmWm& (6.63)
由于
& ( , )
|| & || || || ( , ) ( & , )
( , ) ( & , )( , )
( , ) ( )
( , ) ( )
m W m m W m m
m W m m W m m W m
m W m m m m W m
ij
T
ij ij
T
ij ij
i j S
ij
T
ij ij
T
ij ij
T
ij
i j S
ij
T
ij ij ij ij
T
ij
? +
= + + ? ?
? +
∈
∈
∑
∑
2
2 2 2
2
2
2 2
?
?
(6.64)
由于 &m 与 W 无关,且 ( ) 0, =mm& ,故极小化 (6.64)等价于极小化 (6.65)式
|| || ( , ) ( & , )
( , ) ( )
W m m W m m W mT ij ij T ij ij T ij
i j S
2 2 2? ?
∈
∑
? (6.65)
引入约束条件 Tr[ ]W = 0,利用 Lagrange 乘子法进行极小化,于是目标函数为
J W m m W m m W m WT ij ij T ij ij T ij
i j S
= ? ? +
∈
∑12 22 2|| || ( , ) ( & , ) Tr[ ]
( , ) ( )?
λ (6.66)
将 W 用 W W+ δ 的变分代替,注意到
( ) ( )
( ) ( )( )mWmmWmmWm
mWmWmWmWmW
TTT
TTTTT
dd
ddd
,,2,
,2,
2
2
=
==
(6.67)
从而
( )
( )
???
?
???
?
???
?
???
? +??=
+??=
∑
∑
?∈
?∈
WImmmmmWmmmW
WmWmmmWmmWJ
Sji
T
ij
T
ijijij
T
ij
T
ijij
T
Sji
ij
T
ijijij
T
ijij
T
dl
dldd
)(),(
)(),(
),(Tr
]Tr[,),(
&
&
(6.68)
对于任意的 δW 使上式为 0 的条件为
& ( , )
( , ) ( )
m m W m m m W m m m Iij T T ij ijT ij T ij ij ijT
i j S
? + =
∈
∑
?
λ
(6.69)
对上式两边取迹,注意到 [ ] ( ) [ ] ( )mWmmmWmmmm TTTT ,Tr,,Tr == && ,
[ ] ( ) 1,Tr == mmmmT ,于是有
( & , )
( , ) ( )
m mij ij
i j S∈
∑ =
?
3λ
(6.70)
将 λ 代入式 (6.69)有
[ ] ( )∑∑
?∈?∈
??
?
??
? ?=?
)(),()(),(
,31),(
Sji
T
ij
Sji
T
ijijij
T
ij
T
ijij
T ImmmmmmmWmmmW && (6.71)
将 (6.71)写成分量形式有以下算法。
算法 4.1(利用 N矢量计算速度变换矩阵 )
1. 对给定的光流区域 ?( )S 计算
[ ]
∑ ∑
∑
?∈ =
?∈
??
?
??
? ?=
?=
)(),(
3
1
)(),(
3
1
Syx k
k
xy
k
xyij
j
xy
i
xyij
Syx
l
xy
k
xy
j
xy
i
xy
k
xy
j
xyilijkl
mmmmB
mmmmmmA
&& d
d
(6.72)
其中 i j k l, , , , ,= 1 2 3, ??? ≠
==
ji
ji
ij if0
if1d
, mxy
i
表示 N 矢量 mxy 的第 i 个分量。
2. 速度变换矩阵 W 由下列方程确定
?
?
?
???
?
=
=
∑
∑
=
=
3
1
3
1,
0
i
ii
ij
lk
klijkl
W
BWA
(6.73)
五 . 平面参数与摄象机运动的恢复
注意到速度变换矩阵 W 可以写成对称部分 Ws 和反对称部分 Wa 的和,其中
( ) ( )
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
??+
?=
?
?
?
?
?
?
?
?
?
?
+???
?+?+
?+?
=
1
12
1
100
010
001
3
3)(2)(2)(
2)(3)2(2)(
2)(2)(3)2(
qp
c
b
a
cbaq
p
rr
cqbpa
DAFVEU
FVDACB
EUCBDA
Ws
(6.74)
( ) ( )
?
?
?
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=
?
?
?
?
?
?
?
?
?
?
++
+??
+???
=
0
0
0
121
02)(2)(
2)(02)(
2)(2)(0
12
13
23
ww
ww
ww
qp
c
b
a
cba
r
q
p
r
FVEU
FVCB
EUCB
Wa
(6.75)
算法 4.2 (从速度变换矩阵 W恢复平面参数和摄象机运动 )
1. 设对称矩阵 Ws 的特征值为 σ σ σ1 2 3≥ ≥ ,其所对应的特征向量 u u u1 2 3, , 为相互正
交的单位向量且由于 Tr[ ]W = 0,于是 σ σ σ1 2 3 0+ + = 。
2. 若 σ σ σ1 2 3 0= = = (即 Ws = 0 ),则运动参数为
?
?
?
?
?
?
?
?
?
?
=
?
?
?
?
?
?
?
?
?
?
0
0
0
c
b
a ,
?
?
?
?
?
?
?
?
?
?
?
+?
+
=
?
?
?
?
?
?
?
?
?
?
2)(
2)(
2)(
3
2
1
CB
EU
FV
w
w
w
{ }rqp ,, 不定,否则有以下两组解。
3. 平面的斜率为
p pl q ql= ? ′′ = ? ′′,
其中
332121 uu
l
q
p
ssss ???±=
?
?
?
?
?
?
?
?
?
?
′
′
′
4. r 不定,但与平移速度有如下关系
( )332121 uul
rc
rb
ra
ssss ?+?±′?=
?
?
?
?
?
?
?
?
?
?
5. 旋转速度为
?
?
?
?
?
?
?
?
?
?
×
?
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
?
?
+?
+
=
?
?
?
?
?
?
?
?
?
?
rc
rb
ra
q
p
CB
EU
FV
12
1
2)(
2)(
2)(
3
2
1
w
w
w