摄像机参数标定

数码相机标定

摘要：针对相机标定问题，在做出合理的假设前提下，首先根据相机成像原理，建立了相机的针孔模型，得出了物体的世界坐标与相机的像素坐标之间的转换矩阵M 。这种模型实现简单，但通过已知点的检验，发现这种模型的精度不高。针对针孔线性模型缺点，本文引入了非线性因素（如相机的畸变等），建立了基于Tsai 算法的非线性模型。这种模型需要先用线性模型求出一个较准确的参数，然后根据这些参数利用非线性模型进一步提高精度。针对模型求解过程中涉及到靶标以及靶标像之间的对应关系，确定靶标的像坐标是精确求解模型的关键。针对圆心坐标本文提出了一种求解像坐标圆心的直径交叉法，并用霍夫变换求出圆心相比较，两种算法的相差很小，这种算法理论简单，不需要复杂的数学运算，但精度和霍夫变换相差不多。在问题(3)模型检验中由于图像畸变一般在边缘较大，为了更好体现模型检验的可信性，本文选取图像边缘上的点做检验，计算出模型的最大误差。结果表明非线性模型的精度相对于线性模型精度有很大程度的提高。在问题(4)中，基于以上建立的单相机标定模型的参数，综合考虑两相机标定模型。根据建立的模型中的几何关系，能够得出一般情况下两相机之间的距离和它们的相对角度。当两相机的光学系统光轴之间的夹角近似取零时，即得到特殊的平行视点双相机几何模型，这个模型简单，也适合实际情况。 关键词： 系统标定 霍夫变换 Tsai算法 双目定位

1问题分析

要确定靶标上圆的圆心在像平面的像坐标，即在两者之间建立一种映射关系。容易发现，圆心的像坐标与靶标和相机的相对位置以及相机的内部参数是相关的。因此可以通过设定相关参数及空间相对位置关系得出映射关系。然后利用实际已知的对应点求出某一具体情况下的映射关系，并且进一步利用得到的映射关系求像点坐标。对比求得的像点坐标和实际的像点坐标，确定模型精度。

对于两部相机的标定，可在以上单相机结论的基础上，通过建立与中间点的几何关系，求取两部固定相机的相对位置关系。

2模型假设

理想线性模型的假设：

(1)相机内部参数保持不变； (2)物距远远大于透镜焦距；

(3)物点通过透镜成像后，在像平面上可以找到与之对应的理想像点； 两部相机标定模型的假设：

(1)两部相机型号相同，即内部参数完全相同。 (2)两部相机的光轴共面。

(3)在两部相机拍的图象中，对应靶圆的圆心是匹配的。

3基本符号说明

f —图像平面与光心距离；

dx —每个像素在x 轴方向的物理尺寸；

dy

—每个像素在y 轴方向的物理尺寸；

(u 0, v 0) —相机光轴与图像平面的焦点；

(X w , Y w , Z w )--世界坐标系； (X c , Y c , Z c )--相机坐标系；

R T

--旋转矩阵； --平移向量。

4模型的建立与求解

【1】问题（1）的求解：

数码相机实现了从景物到图像的转换，即实现了从三维景物空间到二维空间中像之间的对应关系，我们主要是求解这种转换的几何关系，也就是三维空间中的点与它在二维空间中的像之间的对应关系。为了定量的描述景物到图像的转换过程，我们首先定义三个坐标系。

相机中的图像是经过处理得到的数字图像，每幅数字图像在计算机中为一个数组，行列的图像中的每个元素即像素是图像点的亮度。在图1-1中，在图像上定义直角坐标系

uO 0v

, 每个像素中的坐标(u , v )分别是该像素在数组中的行数和列数，因此(u , v )是以像

素为单位的图像坐标系的坐标。为了把像素为单位的坐标表示成用物理单位表示出该像素在图像中的位置，需要建立以物理单位（例如毫米）表示的图像坐标系。该坐标系以图像中点O 1为原点，x 轴和y 轴分别平行于u 轴和v 轴。设定在坐标系uO

v

中O 1的坐

标为(u 0, v 0)，每个像素在x 轴与y 轴方向的物理尺寸(每个像素在物理坐标系中的长和宽) 为dx,

图1-1 图像坐标系

u =x /dx +u 0

dy

, 则图像中的任意一个像素在两个坐标系像的坐标的对应关系：

，v

=y /dy +v 0

。

为了分析和计算的方便，把以上变换线性化，用齐次坐标与矩阵的形式表示为：

⎡u ⎤⎡1/dx ⎢⎥⎢v =0⎢⎥⎢⎢⎢⎣1⎥⎦⎣0

01/dy 0

u 0⎤⎡x ⎤

⎥⎢⎥v 0y ⎥⎢⎥1⎥⎦⎢⎣1⎥⎦

其逆关系为：

⎡x

⎤⎡dx

⎥⎢y =0 ⎢⎢⎥⎢⎢⎢⎣1⎥⎦⎣0

0dy 0

-u 0dx ⎤⎡u ⎤

⎥⎢⎥

-v 0dy v

⎥⎢⎥⎥1⎦⎢⎣1⎥⎦

(1)

图1-2相机成像几何关系

相机成像几何关系如图1-2所示，其中O 点为称为相机的光心，X c , Y c 轴与图像的x , y 轴平行，Z c 轴为相机的光轴，它与图像平面的垂直，光轴与图像的交点，即为图像坐标系的原点。由点O c 点与X c , Y c , Z c 轴组成的直角坐标系称为相机坐标系。OO c 为相机光心到像平面的距离。

因为相机可以安放在环境的任何位置，我们需要选择一个基准坐标系来描述相机的位置，并用它描述环境中的任何物体的位置，该坐标系用(X w , Y w , Z w )表示，称为世界坐标系相机坐标系与世界坐标系之间的关系可以用旋转矩阵R 与平移向量T 来描述。因此，世界坐标系和相机坐标系可以用如下关系表示：

⎡X w ⎤⎡X w ⎤

⎡X c ⎤⎢⎥⎢⎥

R T Y Y ⎡⎤⎢⎥⎢w ⎥=M ⎢w ⎥(2) Y c =⎢T

1⎥

⎢⎥

⎢⎣Z c ⎥⎦

⎣0

1⎦⎢Z w ⎥

⎢⎥1⎣⎦

⎢Z w ⎥

⎢⎥1⎣⎦

其中，R 和T 分别是从世界坐标系到相机坐标系的旋转和平移变换，R 为3⨯3正交单位的旋转矩阵，T 为三维平移向量，0T =[0 0 0],M 1为4⨯4矩阵。

为了避免涉及复杂的透镜，需先把描述局限于最简单形式的相机，即针孔相机。在针孔相机中，三维空间中的一个点的图像是从这点发出的并经过针孔的光线与相机后平面的交点，这样就在后平面的底片上产生了图像，从而消除了聚焦的必要性。在像平面上可以找到与之对应的理想像点。

图1-3理想针孔模型

设(X w , Y w , Z w )是三维世界坐标系中目标点P 的三维世界坐标，世界坐标系的原点为O w 。

(X c , Y c , Z c )是同一点P 在相机坐标系中的坐标，相机坐标定义为：中心在O c 点（光学中心），Z c 轴与光轴重合. x , y 是中心在O 点（光轴Z c 与图像平面的交点）平行于X c , Y c 轴的图像坐标系。P 点在图像坐标系与相机坐标系投影的几何关系为：X c

Y c /Z c =y /f

/Z c =x /f

,

, 用齐次坐标表示为：

⎡x ⎤⎡f

⎢⎥⎢y =0⎢⎥⎢⎢⎣1⎥⎦⎢⎣0

0f 0

001

⎡X c ⎤

0⎤⎢⎥⎥⎢Y c ⎥0

⎥⎢Z ⎥

c

0⎥⎥⎦⎢

1⎣⎦

Z c

(3)

我们得到四个坐标系之间的三个关系：公式（1）表示了图像坐标系与像素坐标系之

间的矩阵关系；公式（2）表示了世界坐标系与相机坐标系之间的矩阵关系; 公式（3）表示了相机坐标系与图像坐标系之间的矩阵关系。

为了得到像素坐标系与世界坐标系之间的矩阵关系，将公式（1），（2）式代入（3）式可以得到（4）式，即世界坐标系与像素坐标系之间的矩阵关系：

⎡X w ⎤

u a 0u 0⎡⎤⎡x ⎤0⎢⎥ T ⎤Y w ⎢⎥⎢⎥⎡R ⎢⎥Z c v =0a y v 00⎢T =M 1M 2X =MX ⎥ ⎢⎥⎢⎥0⎢⎥1Z ⎣⎦w ⎢⎢010⎥ ⎢⎥⎣1⎥⎦⎣0⎦

1⎣⎦

将上式写成：

⎡u ⎤⎡m 11⎢⎥⎢Z c v =m 21

⎢⎥⎢⎢⎣1⎥⎦⎢⎣m 31

m 12m 22m 32

m 13m 23m 33

⎡X w ⎤

m 14⎤⎢⎥

⎥⎢Y w ⎥m 24

⎥⎢Z ⎥

w

m 34⎥⎥⎦⎢

⎣1⎦

(4)

式中，(X w ,Y w ,Z w , 1)为空间三维点的世界坐标，(u , v , 1)为图像坐标，m ij 为透视矩阵M 的 元素，包含三个方程，为了消除

m 11X m 21X

w

Z c

，整理方程后可得关于

w

m ij

的线性方程。

+m 12Y w +m 13Z w +m 14-uX

w

m 31-uY w m 32-uZ m 33=um

34

(5) (6)

w

+m 22Y w +m 23Z w +m 24-vX

w

m 31-vY w m 32-vZ w m 33=vm 34

（5）（6）式描述了世界坐标系中的点与相应的像素坐标系的点的对应关系。只要知道坐标系中的一个点的坐标值和与之对应于像素坐标系中的坐标值，就可以得到（5）（6）两个方程。取n 个点的世界坐标值和像素坐标值，就可以得到2n 个方程，下面用矩阵的形式表示这些方程：

⎡m 11⎤⎡u 1m 34⎤⎢⎥⎢⎥m 12v 1m 34⎢⎥⎢⎥⎢m 13⎥⎢⎥⎢⎥⎢⎥

-u 1Z w 1⎤⎢m 14⎥⎢⎥

⎥⎢⎢ ⎥m 21⎥-v 1Z w 1⎥

⎢⎥⎢⎥⎥⎢m 22⎥ ⎢⎥

=⎥⨯⎢⎥⎢ ⎥

⎥⎢m 23⎥⎢⎥⎥-u n Y wn ⎢m 24⎥⎢ ⎥⎥⎢⎥⎢⎥

-v n Z wn ⎥⎦⎢m 31⎥⎢⎥

⎢m 32⎥⎢⎥⎢⎥⎢⎥

u m m ⎢33⎥⎢n 34⎥

⎢m ⎥⎢v n m 34⎥

⎣⎦⎣34⎦

⎡X w 1

⎢⎢0⎢⎢⎢⎢X

wn ⎢⎢⎣0

Y w 20

Z w 30... ...

10... ... 10

0X

w 1

0Y w 1

0Z w 1

01

-u 1X -v 1X

w 1w 1

-u 1Y w 11-v 1Y w 1

... ... 0X

wn

(7)

Y wn 0

Z wn 0

0Y wn

0Z wn

01

-u n X -v n X

wn wn

-u n Y wn -v n Y wn

M

矩阵乘以任意不为零的常数并不影响(X w , Y w , Z w ) 与(u , v )的关系，因此，在上式中可

=1，从而得到关于M

以指定m 34

矩阵其他元素的2n 个线性方程，这些未知元素的个数

KM =U

为11个，记为11维向量M ，上式简写成：

(8) 其中，K 为上式(7)左边2n ⨯11矩阵，m 为未知的11维向量，U 为上式(7)右边的2n 维向量；K ，U 为已知向量。当2n >11时，我们可用最小二乘法求出上述线性方程的解：

-1T T (9) M =(K K )K U

从而求得投影矩阵M （ Matlab程序在附录四）。对于标靶上圆的圆心在该相机平面上的像坐标，就可以通过投影矩阵M 将标靶上圆心的世界坐标值映射到平面上的像素坐标，进而得到圆心的像素坐标值。 【2】问题（2）的解答

（1）特征值的提取（直径交叉法）

相机图像可以看作像素的集合体，我们可以对像素进行运算。将图像看作一个二维矩阵，矩阵中的元素为像素，通过选择矩阵中的元素的值可以检测出图像的某种特性。对于标靶的像中5个圆称为孤立圆，孤立圆的强度与背景像素不同，因此可以检测出这

些孤立圆的圆心的像素。5个圆的圆心坐标并不能求解投影矩阵M ，需要根据两幅图像中的比例关系再提取一个点，我们认为标靶示意图中对角线的交点P6（中心点）与标靶像中P1P5和P3P4的交点P6对应的可能性最大，误差最小。因此，取P6为个圆心点。靶标图像中6个点的位置如图2-1所示：

图2-1靶标图像中的6个像点

要确定圆心的像素坐标（Matlab 程序在附录二）：

a. 解决图像边界的问题，题目中给出了相机分辨率为1024*768（4/3的比例关系），设定在像素坐标系中像素的有效范围:u

∈(0, 1024)

, v ∈(0, 768) 。

b. 对图像进行预处理。将标靶的像通过计算机传入Matlab 中，每一个孤立圆的边界用相隔均匀的四个点来标定，四个点的像素坐标通过Matlab 的函数（impixel ）可得，从而可以初步缩小圆心的提取范围。

c. 确定圆心的像素坐标。由imread 函数将图像中的像素读出，为了明显观察像素的灰度值差异，设定背景像素灰度值为0，圆的像素灰度值为1。在圆心的提取范围内，1的数目最多的行和列的交点即为圆心的位置，从而可以确定圆心的像素坐标。

考虑到霍夫变换是图像处理中从图像中识别几何形状的基本方法之一，应用广泛。基本思想就是把图像平面上的点对应到参数平面上的线，最后通过统计特性来解决问题。假如图像平面上有两条直线，那么最终在参数平面上就会看到两个峰值点，最基本的霍夫变换是从黑白图像中检测直线(线段) 。针对本题中的靶标图像圆心的确定，即从一幅图像中检测出半径以知的圆形来。我们可以取和图像平面一样的参数平面，以图像上每一个前景点为圆心，以已知的半径在参数平面上画圆，并把结果进行累加。最后找出参数平面上的峰值点，这个位置就对应了图像上的圆心(如图2-2所示) 。在这个问题里，图像平面上的每一点对应到参数平面上的一个圆。

图2-2 霍夫算法提取的圆心

通过Matlab 编程，程序（附录一）我们得到了圆心的像素坐标值如下：

表二：标靶像中5个圆心的像素坐标

下面我们对上述两种方法加以比较。

交叉直径法利用图像全局特性而将像素集合为一个特定区域，没有考虑到相机在成像的过程中，由于内部参数和外界噪音的因素而造成图形的畸变。霍夫变换对噪声和线间断的影响较小。为了减少噪音等因素对图像造成的影响，我们采用霍夫变换求得的结果。

现在我们得到了6个圆心的世界坐标值和对应的像素坐标值，每一圆心对应两个方程式(5)(6),用最小二乘法求解方程式(9)可得：

⎡3. 1064

⎢

M =0. 0612

⎢⎢⎣-0. 001

0. 3970-3. 1706-0. 001

231. 0096175. 52340. 0138

231. 0096⎤

⎥

175. 5234

⎥

1. 0000⎥⎦

(10)

【3】问题（3）的解答

基于霍夫变换是将不连续的边缘像素点连接起来，组成一个封闭的边界，边缘点的精度和稳定性是我们关心的重点, 图像畸变在边缘上一般较大。对模型的检验，我们采取边缘检测的方法。取P1和P3的中点P7，P1和P4的中点P8，P4和P5的中点P9，P3和P5的中点P10以及P1、P2、P3、P4、P5共9个边缘检验点。标靶像中边缘点的坐标如下图3-1所示：

图3-1靶中边缘点的坐标

通过投影矩阵M 可得到了转换之后10个边缘点的像素坐标。我们设定P7 、P8、P 9的实际的像素坐标是通过表二中前5个点的像素坐标经过适当的比例运算得到的。每个点转换后的像素坐标(

u , v

'

'

)和实际像素坐标(u , v )之间的误差用相对误差δ（Matlab

程序在附录三）来表示：

δ=将(

u , v

'

'

((

u -u

'

)2

+v -v

(

'

))/(u

2

2

+v

2

) (11)

)、(u , v )、δ的求解结果如下：

误差的分析与模型的改进：

由上表可以看出部分P 点的坐标误差相对较大，原因经分析知：相机的几何模型，只是一种理想的线性模型。当用实际的相机拍摄景物的图像时，由于镜头系数的精密度有限，电子扫描存在非线性、光敏元件矩阵分布不可能绝对一致，因此不可避免的存在几何畸变。在前面的理想线性模型（针孔模型），相机标定除了考虑主点坐标，焦距和歪斜系数S 等线性畸变，还要考虑透镜畸变也会造成成像坐标的偏差，包括透镜的径向畸变和切向畸变。 非线性畸变的分析:

非线性畸变包括镜头的径向畸变和切向畸变, 两种畸变比较, 径向畸变为影响精度的主要原因. 所以, 主要考虑径向透镜畸变, 并与一个二阶多项式近似表示(径向畸变为一个高阶多项式):

x d =x u 1+kr d

(

2

), y

d

=y u 1+kr d

(

2

)

式中, r d

2

=x d +y d ;

22

k 为畸变系数。

针对针孔透镜模型没有考虑畸变等非线性因素所造成的偏差, 我们提出了考虑径向畸变的Tsai 非线性模型，这种模型的求解需要给定设定的初值，如果给定的初值不准确将很难得到正确的结果。因此我们考虑将前面的针孔线性模型与Tsai 非线性模型相结合，先用线性模型求解相机的参数，再以解的参数作为初始值，考虑畸变因素，利用非线性优化方法，进一步提高标定精度，具体算法的实现步骤如下：

选取世界坐标系使：Z w =0。

第一步：求解旋转矩阵R 和T 中的t x , t y ；

（1） 确定n 个图像点的世界坐标系和图像坐标系；

（2） 对每个点，联立N 个方程，利用最小二乘法求解超定方程，可求得以下变量：

r 1=r 1/t y

'

，r 2' =r 2/t y ，t x =t x /t y

' ，r 4' =r 4/t y ，r 5=r 5/t y

'

。

利用R 的正交性可以算出t y ，和r 1--r 9； 第二步：计算有效焦距

f

, T 的t x 分量和畸变系数k ；

x u =f ⨯X /Z =x d 1+kr d

(

2

)；

2

y u =f ⨯Y /Z =y d 1+kr d 。

()

式中，r d 2设H x

=x d +y d

22

，待求变量

，H

H

x

y

f

, k , t z 。

，W

=r 7x w +r 8y w

=r 1x w +r 2y w +t x =r 4x w +r 5y w +t y

2

，

f k =f ⨯k

，可得：

⨯f +H

x

⨯r d ⨯f k -y d ⨯t z =y d ⨯W

2

H

y

⨯f +H

y

⨯r d ⨯f k -y d ⨯t z =y d ⨯W

。

对N 个特征点，利用最小二乘法对上述两个方程进行联合最有参数估计，就可得，进而求得f ，f k ，t z 。

根据附录五的Tsai 算法的Matlab 程序算得：

⎡0. 9484⎢

R =0. 0253

⎢⎢⎣0. 3162

-0. 3109-0. 12330. 9424

⎤

⎥⎥⎥⎦

⎤

⎥

-0. 9921

⎥

-0. 1090⎥⎦0. 0628

；

⎡7. 9924

⎢

T =-35. 0115

⎢⎢⎣1. 4545

f =45. 8957

；

；

-5

k =-3. 9012⨯10

。

从而可以确定出相机的内外参数N 。

【4】问题（4）的解答

两部相机标定模型

问题中要求用靶标给出两部固定相机的相对位置，即已知靶标上的目标点的世界坐标和用两部固定相机中拍摄出的像坐标来确定两部固定相机的相对位置。如图4-1所示，D 在两个CCD 图像上的成像位置不同，进而根据这个差异和D 点在世界坐标系中的坐标来确定两个固定相机之间相对位置。

图4-1 双相机光轴共面定位模型图象

图4-1中，CCD1、CCD2为两个位置相对固定的数码相机的像平面；

v 1、v 2

分别为两个光学系统的像距；

d 1、d 2分别为两个光学中心到所要求的对称中心的距离； L 1、L 2分别为两个光学系统的物距；

θ1、θ2分别为两个光学系统光轴与所要求的对称线的夹角。 空间坐标系中的点在O 1X 1Y 1Z 1坐标系下的几何光学投影为:

⎡x 1⎤

⎡1

00

⎤⎡x ⎤⎡0⎤

⎢-y 1⎢⎢y ⎥1⎥=d 1)1sin θ-y sin θ1+z cos θ1⎢0cos θ1sin θ⎥⎢⎥⎢⎥1⎥⎢y -d 1⎥-⎢0

⎥ ⎢⎣z 1⎥⎦

⎢⎣0

-sin θ1

cos θ1⎥⎦⎢⎣z ⎥⎦⎢⎣-v 1⎥⎦在O 2X 2Y 2Z 2坐标系下的几何光学投影为:

⎡x 2⎤

⎡1

00

⎤⎡x ⎤⎡0⎢⎥

-y 2⎢y 2⎥=d -y sin θ2+z cos θ2

)⎢2sin θ2

⎢0cos θ2sin θ⎥⎤⎥⎢⎢y -d ⎥⎢⎥22⎥-⎢0

⎥

⎢⎣z 2⎥⎦

⎢⎣

0-sin θ2

cos θ⎥⎢2⎦⎣z ⎥⎦⎢⎣-v 2⎥⎦由(15)(16)两式可得:

d 1sin θ1x 1-y sin θ1x 1+z cos θ1x 1+v 1x =0 (14)

12）（

(13)

d 2sin θ2x 2-y sin θ2x 2+z cos θ2x 2+v 2x =0 (15)

P -P 2y =0 （16） 1P 3+P 2z =0 (17) d =d 1+d 2 (18)

θ=θ1+θ

2

（19）

其中:

P 1=(y 1d 1tan θ1-v 1d 1)(y 2-v 2tan θ2)-(y 1+v 1tan θ1)(y 2d 2tan θ2+v 2d 2)

(20) (21)

P 2=(y 1tan θ1-v 1)(y 2-v 2tan θ2)-(y 1+v 1tan θ1)(y 2d 2tan θ2+v 2d 2)

P 3=(y 1tan θ1-v 1)(y 2d 2tan θ2+v 2d 2)-(y 1d 1tan θ1-v 1d 1)(y 2d 2tan θ2+v 2) （22）

给定P 点在三个坐标系中的坐标(x , y , z )、(x 1, y 1, z 1)、(x 2, y 2z 2), 由式（14）-（22）就可解得θ1, θ2, d 1, d 2的值, 进而得出d 和θ的值，即两部固定相机的相对距离和相对角度。特别地，当θ无限趋向于零时，即得到平行视点双目立体几何模型。

5模型的评价

针孔模型虽然做了理想的简化假设，但这个模型对大多数通用透镜来说也相当适合，最重要的是，它简单实用而不失准确性，其误差在合理的范围之内。然而，对于高精度要求的系统标定来说，针孔模型达不到要求。因此，在针孔线性模型的基础上，通过引入歪斜系数建立的线性模型在一定程度上提高了模型精度。为了使模型更贴近实际情况，又建立了非线性的畸变模型。本文建立Tsai 的畸变模型，考虑了畸变中的主要因素，避免引入过多的非线性系数，既提高了精度，又具有较好的实时性。

尽管前面所提出的各种标定方法在理论上都是切实可行的，但仍有许多地方需要做进一步的研究。

（1）图象一般包含有噪声，即使这种噪声非常小，内部参数的实际解与由约束关系所得解之间的差异仍然是相当大的，怎样提高解的鲁棒性成为相机标定领域的一个必须解决的问题。

（2）相机标定所解决的问题归根揭底是得到一组非线性方程的解，一般情况下使用各种优化方法，但实际上所得优化解往往不是全局的，因此，进一步研究更合理的求解非线性方程的方法同等重要。

（3）不确定性是相当重要的，因为它决定计算参数的可信度，相机定标参数的不

确定性对三维重建有何影响和约束关系的不确定性是怎样传播等这些都有待进一步深入研究。

参考文献

[1] 罗均，谢少荣等编著，特种机器人，北京：化学工业出版社，2006 [2] 方建军，何广平编著，智能机器人，北京：化学工业出版社，2003

[3] [德]马科斯. 玻恩，[美]埃米尔. 沃尔夫编著，杨葭蓀等译，光学原理，

北京：电子工业出版社，2005

[4] 廖延彪编著，光学原理与应用，北京：电子工业出版社，2006

[5] 张可，张锦好，双目立体视觉系统的非线形模型建立，传感器与微系统：

2006第25卷 第10期 [6] 郝晓辰等，双目立体视觉测量模型与同名点匹配研究，仪器仪表学报：2005

第26卷 第8期增刊

附录

附录一: 霍夫变换提取圆心； c=imread('c.bmp' ); a=rgb2gray(c); tic;

[d,e,f] = CircularHough_Grd(a, [20 30],5,50); toc;

figure(1); imagesc(a); colormap('gray' ); axis image ; title(' 提取圆心' ); hold on ;

plot(e(:,1), e(:,2), 'r+'); for k = 1 : size(e, 1),

DrawCircle(e(k,1), e(k,2), f(k), 32, 'b-' ); end

附录二： 直径交叉法提取圆心；

function f=tiqu1(m1,m2,n1,n2) f=zeros(1,2);

a=imread('c.bmp' ); b=zeros(1,m2-m1+1); c=zeros(1,n2-n1+1); for i=m1:m2

for j=n1:n2

if (a(i,j)

b(1,i-m1+1)=b(1,i-m1+1)+1; end end end

for i=m1:m2

for j=n1:n2

if (a(i,j)

c(1,j-n1+1)=c(1,j-n1+1)+1; end end end [e,k]=max(b);

[g,h]=max(c); f(1,1)=k+m1; f(1,2)=h+n1; f

附录三：

求取数据的相对误差； x1=[303.9711

481.4348,611.4297,433.9660,422.7229,323.0290,639.8407,583.0187,284.9132]; x2=[306.8180,481.3080,617.2680,442.778,419.7412,326.5492,637.1892,597.4892,286.8492];

y1=[345.667 201.4409 358.2630 502.4588 197.0116 189.4843 213.3976 501.7492 503.1684];

y2=[348.8382,192.5482 354.6882 510.9782 191.2928 189.4568 195.5766 512.6368 506.5168]; n=size(x1,2); f1=0;

f2=zeros(1,n); for i=1:n

f1=((x1(i)-x2(i)).^2+(y1(i)-y2(i)).^2).^0.5; f2(i)=f1/(x1(i).^2+y1(i).^2).^0.5 end

f=sum(f2); f=f/n; 附录四：

M 矩阵的求解;

k=[-50 -50 1 1 0 0 0 0 319*50 319*50 -319*1; 0 0 0 0 -50 -50 1 1 50*184 50*184 -184*1; -20 -50 1 1 0 0 0 0 20*425 50*425 -425*1; 0 0 0 0 -20 -50 1 1 193*20 193*50 -193*1; 50 -50 1 1 0 0 0 0 -50*635 50*635 -635*1; 0 0 0 0 50 -50 1 1 -212*50 50*212 -212*1; -50 50 1 1 0 0 0 0 280*50 -50*280 -280*1; 0 0 0 0 -50 50 1 1 502*50 -50*502 -502*1; 50 50 1 1 0 0 0 0 -582*50 -582*500 -582*1; 0 0 0 0 50 50 1 1 -503*50 -503*50 -503*1; 0 0 1 1 0 0 0 0 0 0 -535*1; 0 0 0 0 0 0 1 1 0 0 -239*1;];

u=[319 184 425 193 635 212 280 502 582 503 432 344]; m=pinv(k)*u';

m1=[m(1,1:4);m(1,5:8);m(1,9:11),1]; 附录五：

利用Tsai 算法求解空间转换矩阵R 和位移矩阵T 以及焦距

f

和畸变系数k ；

clc clear m1=[-50 50 -20 50 50 50 -50 -50 50 -50 ];

m2=[322.8870 189.3366 422.7260 196.9228 639.8416 213.3503 582.7751 502.9878 284.923 501.7052 ];

xw=m1(:,1); yw=m1(:,2); Xf=m2(:,1); Yf=m2(:,2); [M,N]=size(xw); zw=zeros(M,1); Ncx=1024; Nfx=768; Cx=1024/2; Cy=768/2; dx=1/3.78; dy=1/3.78; sx=1;

[R, T, f, k1] = Tsai(Xf, Yf, xw, yw, zw, Ncx, Nfx, dx, dy, Cx, Cy, sx) % R 旋转矩阵,T 平移向量 f 是焦距长度, k1是畸变系数 %在同一个平面上,ZW 相同, 都取10

function [R, T, f, k1] = Tsai(Xf, Yf, xw, yw, zw, Ncx, Nfx, dx, dy, Cx, Cy, sx) % Stage 1 --- Compute 3D Orientation, Position (x and y): % a) Compute the distored image coordinates (Xd, Yd) Procedure: dxp = dx * Ncx / Nfx;

X = Xf - Cx; Y = Yf - Cy;

Xd=sx^(-1)*dxp*(Xf-Cx); Yd=dy*(Yf-Cy);

% b) Compute the five unknowns Ty^(-1)*r1, Ty^(-1)*r2, Ty^(-1)*Tx, Ty^(-1)*r4, Ty^(-1)*r5

% r1p=Ty^(-1)*r1; % r2p=Ty^(-1)*r2; % Txp=Ty^(-1)*Tx; % r4p=Ty^(-1)*r4; % r5p=Ty^(-1)*r5;

A=[Yd.*xw Yd.*yw Yd -Xd.*xw -Xd.*yw]; B=Xd; C=A\B; r1p=C(1); r2p=C(2); Txp=C(3); r4p=C(4); r5p=C(5); clear A B C ;

% c) Compute (r1,...,r9,Tx,Ty) from (Ty^(-1)*r1, Ty^(-1)*r2, Ty^(-1)*Tx, Ty^(-1)*r4, Ty^(-1)*r5):

% 1) Compute |Ty| from (Ty^(-1)*r1, Ty^(-1)*r2, Ty^(-1)*Tx, Ty^(-1)*r4, Ty^(-1)*r5): C=[r1p, r2p; r4p, r5p];

Sr=r1p^2 + r2p^2 + r4p^2 + r5p^2; if rank(C)==2

Ty2=( Sr - (Sr^2-4*(r1p*r5p-r4p*r2p)^2)^(0.5) )/(2*(r1p*r5p-r4p*r2p)^2); else

z = C(abs(C) > 0);

Ty2 = 1.0 / (z(1)^2 + z(2)^2); end

Ty = sqrt(Ty2); clear C Sr Ty2 z

% 2) Determine the sign of Ty: [ymax i] = max(Xd.^2 + Yd.^2); r1 = r1p*Ty; r2 = r2p*Ty; r4 = r4p*Ty; r5 = r5p*Ty; Tx = Txp*Ty;

x = r1*xw(i) + r2*yw(i) + Tx; y = r4*xw(i) + r5*yw(i) + Ty;

if (sign(x) == sign(Xf(i))) & (sign(y) == sign(Yf(i))), Ty = Ty; else

Ty = -Ty; end

clear ymax i x y

% 3) Compute the 3D rotation matrix R, or r1, r2,...,r9 r1 = r1p*Ty;

r2 = r2p*Ty; r4 = r4p*Ty; r5 = r5p*Ty; Tx = Txp*Ty;

s = -sign(r1*r4 + r2*r5);

R=[r1, r2, (1-r1^2-r2^2)^(0.5); r4, r5, s*(1-r4^2-r5^2)^(0.5)]; R = [R(1:2,:); cross(R(1,:), R(2,:))];

r7 = R(3,1); r8 = R(3,2); r9 = R(3,3);

y = r4*xw+r5*yw+Ty; w = r7*xw+r8*yw;

z = [y -dy*Y] \ [dy*(w.*Y)]; f = z(1);

if f

R(1,3) = -R(1,3); R(2,3) = -R(2,3); R(3,1) = -R(3,1); R(3,2) = -R(3,2); end

r3 = R(1,3); r6 = R(2,3); r7 = R(3,1); r8 = R(3,2); clear s y w z

% 2) Stage 2 --- Compute Effective Focal Length, Distortion Coefficients, and z Position: % d) Compute an approximation of f and Tz by ignoring lens distortion: y = r4*xw+r5*yw+Ty; w = r7*xw+r8*yw;

z = [y -dy*Y] \ [dy*(w.*Y)]; f = z(1); Tz = z(2);

% Compute the exactly solution for f, Tz, k1: params_const = [r4 r5 r6 r7 r8 r9 dx dy sx Ty]; params = [f, Tz, 0]; % add initial guess for k1

[x,fval,exitflag,output] = fminsearch( @Tsai_8b, params, [], params_const, xw, yw, zw, X, Y);

f = x(1); Tz = x(2); k1 = x(3);

T=[Tx, Ty, Tz]';

% fval the value of the objective function fun at the solution x. fval

% exitflag that describes the exit condition of fminsearch % >0 Indicates that the function converged to a solution x.

% 0 Indicates that the maximum number of function evaluations was exceeded. %

% output that contains information about the optimization % output.algorithmThe algorithm used

% output.funcCountThe number of function evaluations % output.iterationsThe number of iterations taken output

function f = Tsai_8b(params, params_const, xw, yw, zw, X, Y) % unpack the params f = params(1); Tz = params(2); k1 = params(3); % unpack the params_const r4 = params_const(1); r5 = params_const(2); r6 = params_const(3); r7 = params_const(4); r8 = params_const(5); r9 = params_const(6); dx = params_const(7); dy = params_const(8); sx = params_const(9); Ty = params_const(10);

rsq = (dx*X).^2 + (dy*Y).^2;

res = (dy*Y).*(1+k1*rsq).*(r7*xw+r8*yw+r9*zw+Tz) - f*(r4*xw+r5*yw+r6*zw+Ty); f = norm(res, 2);