张氏标定法原理及其改进1

张正友算法原理及其改进

由于世界坐标系的位置可以任意选取，我们可以假定世界坐标系和摄像机坐标系重合，故定义模板平面落在世界坐标系的0W Z =平面上。用i r 表示R 的每一列向量，那么对平面上的每一点，有：

[][]123120111W W W W X u X Y s v A r r r t A r r t Y ⎡⎤

⎡⎤

⎡⎤⎢⎥

⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥

⎢⎥⎢⎥⎢⎥⎣⎦

⎣⎦

⎣⎦

(3.1) 这样，在模板平面上的点和它的像点之间建立了一个单应性映射H ，又称

单应性矩阵或投影矩阵。如果已知模板点的空间坐标和图像坐标，那么就已知m 和M ，可以求解单应性矩阵H 。)1,,(w w Y X )1,,(v u

因为11W W u X s v H Y ⎡⎤⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦,其中11

121321222331

32

1h h h H h h h h h ⎡⎤

⎢⎥=⎢⎥⎢⎥⎣⎦

,可推出： 111213

21222331321

W W W W W W su h X h Y h sv h X h Y h s h X h Y =++⎧⎪

=++⎨⎪=++⎩

(3.2) 故，

1112133132212223313211W W W W W W W W h X h Y h u h X h Y h X h Y h v h X h Y ++⎧=⎪++⎪

⎨

++⎪=⎪++⎩

(3.3) 将分母乘到等式左边，即有

3132111213

3132212223W W W W W W W W uX h uY h u h X h Y h vX h vY h v h X h Y h ++=++⎧⎨++=++⎩ (3.4)

又令[]T

h h h h h h h h h 323123

2221131211

='，则

1

0000

01W W W W W

W

W

W X Y uX uY u h X Y vX vY v --⎡⎤⎡⎤

'=⎢⎥⎢⎥--⎣⎦⎣

⎦

(3.5) 多个对应点的方程叠加起来可以看成Sh d '=。利用最小二乘法求解该方程，即1()T T h S S S d -'=，进而得到H 。 摄像机内部参数求解

在求取单应性矩阵后，我们进一步要求得摄像机的内参数。首先令i h 表示H 的每一列向量，需要注意到上述方法求得的H 和真正的单应性矩阵之间可能相差一个比例因子，则H 可写成：

[][]12312h h h A r r t λ= (3.6)

又因为1r 和2r 是单位正交向量，所以有

1210T T h A A h --= (3.7)

11

1122T T T T h A A h h A A h ----= (3.8)

这样就为内参数的求解提供了两个约束方程。 下面，令

()()()111213

1212223

313233

002222000222222

222

000000002222222111T B B B B A A B B B B B B v u v u v v u v u v u v v γβγααβαβγγβγγαβ

αββαββγγβγβγβαβαββαββ--⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦⎡⎤--⎢⎥

⎢⎥⎢⎥--⎢⎥=+--⎢⎥⎢⎥---⎢⎥--++⎢⎥⎣⎦

注意到B 是一个对称矩阵，所以它可以由一个6维向量来定义，即

[]11

12

22

13

23

33T

b B B B B B B = (3.9)

设H 的第i 列向量为123[,,]T i i i i h h h h =，因此有

T T i j ij h Bh V b = (3.10)

其

中

1

,

,

T

i j V h ⎡⎤=+

+

+

⎣

⎦，那么，

就可以将内参数的两个约束写成关于b 的两个方程为：

1211220T

T T V b V V ⎡⎤

=⎢⎥

-⎢⎥⎣⎦

(3.11) 如果有n 幅图像的话，把它们的方程式叠加起来，得到

0Vb = (3.12)

其中，V 是一个26n ⨯的矩阵。当3n ≥时，一般情况下，b 可以在相差一个尺度因子的意义下唯一确定；当2n =时，此时的方程的个数少于未知数的个数，我们可以加上一个附加约束0γ=，即120B =，因此可用[]0100000b =作为式(3.12)的一个附加方程。方程(3.12)的最小二乘解即是T V V 的最小特征值对应的特征向量，将该向量归一化即得到要求的b ，进而得到B ；当1n =时，两个方程只能解两个未知数，我们可以假定光心投影在图像的中心，从而求出摄像机在水平和垂直方向上的最大倍数。

一旦b 被求出，就能根据下面两种方法计算出摄像机的内参数矩阵A ： (1)由b 构造出B ，再利用Cholesky 矩阵分解算法求解出1A -，再求逆得到A 。

(2) 由b 构造出B ，在相差一个尺度因子的意义下(1T B A A λ--=,其中为λ尺度因子) ,由绝对二次曲线的性质，很容易求出摄像机的内部参数：

2012131123112212()/()v B B B B B B B =-- 2221301213112311[()]/B B v B B B B B λ=-+-

x f =

y f =212/x y B f f αλ=-

20013//x x u v f B f αλ=- 摄像机外部参数求解

由每幅图像的单应性矩阵H 和上一节的计算结果就可以求得每幅图像的外部参数。

一旦A 求得后，根据式(3.6)，每幅图像的外部参数很容易求出：

111r A h λ-= 122r A h λ-= 312r r r =⨯ 13t A h λ-=

这里的尺度因子11121/1/A h A h λ--==。当然，由于图像必然有噪声，因此这样解得的123[,,]R r r r =并不能完全满足旋转矩阵的性质，所以要从一个给定的矩阵求解一个最佳的旋转矩阵。 非线性优化（优化内参）

以上我们所得到的摄像机的内参数矩阵和每幅图像对应的外参数矩阵都只是一个粗糙解，没有具体的物理意义，可以通过最大似然估计对所有参数进行非线性优化，进一步求精。在这里可以假定有n 幅关于模板平面的图像，模板平面上有m 个标定点，那么可建立评价函数：

()2

11

,,,n

m

ij i j i j C m m A R t M ===-∑∑ (3.13)

其中ij m 是第i 幅图像中的第j 个像点，i R 是第i 幅图坐标系的旋转矩阵，i t 是第i 幅图坐标系的平移向量，j M 是第j 个点的空间坐标，(),,,i j m A R t M 是通过这些已知量求得的像点坐标。

由于旋转矩阵有9个参量但是只有三个自由度，因此可用三个参量的矢量来表示，即一个旋转可由一个三维向量即旋转向量来表示，他的方向就是旋转轴的方向，他的模等于旋转角。

由三个欧拉角参数确定，123(,,)T r r r r 是旋转矩阵的罗德里克(Rodrigues)表示，R 与r 之间的关系由公式给出：

22sin /[](1cos )/[]x x R I r r θθθθ=++- (3.14)

其中，旋转向量T z y x r ),,(=，定义由它构成的反对称矩阵为