张正友相机标定算法学习笔记
张正友标定算法学习笔记
常见的相机标定方法包括基于控制场的标定方法、基于自检校模型的光束法平差、及基于二维平面格网的标定方法等。
自标定方法在不需建立标定场的情况下,利用目标周围的影像与影像之间的对应关系对相机进行标定,该方法较为灵活且成本低,缺点是稳健性较差。比较有代表性的是张正友标定法。
基本流程
1. 获取观测数据
对棋盘格拍照,则得到多张图像的多个角点的像点和世界坐标系坐标对 m ^ \hat m m^和 M ^ \hat M M^
M
^
=
[
X
Y
Z
1
]
\hat M=\begin{bmatrix} X\\ Y\\ Z\\ 1 \end{bmatrix}
M^=⎣⎢⎢⎡XYZ1⎦⎥⎥⎤
其中,
Z
=
0
Z=0
Z=0,因为棋盘格是平面的。
m
^
=
[
μ
ν
1
]
\hat m=\begin{bmatrix} \mu\\ \nu\\ 1 \end{bmatrix}
m^=⎣⎡μν1⎦⎤
m ^ = 1 z c A [ R t ] M ^ (1) \hat m=\frac{1}{z_c} A[R \quad t]\hat M \tag{1} m^=zc1A[Rt]M^(1)
拍照时候的注意事项和疑问:
- 标定板的平整度,棋盘格大小(尽量大),图像画质,角点分布是否均匀,角点覆盖整个画面。
- 棋盘格固定,改变相机的姿态和距离?距离和角度范围怎么确定呢?
- 要得到标定板角点的物理坐标值,需要根据已知的棋盘格大小和世界坐标系原点。这个怎么设置呢?
2. 计算单应性矩阵
计算 m ^ \hat m m^与 M ^ \hat M M^对应的单应性矩阵 H H H
z
c
m
^
=
H
M
^
(2)
{z_c} \hat m = H \hat M \tag{2}
zcm^=HM^(2)
其中,
H
H
H包含了相机内参矩阵
A
A
A、外参旋转矩阵
R
R
R、外参平移参数
t
t
t,
H
=
A
[
R
t
]
(3)
H = A[R \quad t] \tag{3}
H=A[Rt](3)
A
=
[
f
x
γ
μ
0
0
f
y
ν
0
0
0
1
]
A = \begin{bmatrix} f_x& \gamma & \mu_0\\ 0& f_y & \nu_0 \\ 0& 0& 1 \end{bmatrix}
A=⎣⎡fx00γfy0μ0ν01⎦⎤
其中,
f
x
,
f
y
f_x, f_y
fx,fy分别表示两个方向的焦距,
μ
0
,
ν
0
\mu_0, \nu_0
μ0,ν0表示光心点的坐标,即分别表示相机感光板中心在像素坐标系下的坐标,
γ
\gamma
γ表示两个方向之间的倾斜系数,一般来说默认等于0。
注意这个步骤只是解算了 H H H, H H H是齐次矩阵,有8个独立未知元素。当一张图片上的标定板角点数量等于4时,即可列方程按照最小二乘求解该图片对应的矩阵 H H H。
另外一个说法: H H H可以根据棋盘格角点的实际坐标和图像对应坐标求出,参考多目几何直接线性化方法(DLT),每一张图像都会形成一个单应矩阵。
3. 由单应性矩阵计算内参
对 H H H分解,计算 B = A − T A − 1 B=A^{-T}A^{-1} B=A−TA−1矩阵,并利用 B B B求解 A A A内参矩阵
这里就是一些关于矩阵的公式推导,要点有:
- 旋转矩阵的特点(r1, r2 是单位正交的向量)。
- B 是一个对称矩阵,有6个参数需要求解,一个棋盘格有两个方程,所以至少需要拍摄3张不同位置或者角度棋盘格图像(事实上一般需要15到20张标定板图片),列方程按照最小二乘求解。
- 需要需要拍摄10几张以上的像片,利用坐标点对。
- 利用 B B B求解 A A A内参矩阵的元素是直接推导出来的。
- 具体的公式推导参考这个博客.
4. 由内参推导各图像的外参
根据 A A A计算,每张棋盘格的外参 R R R和 t t t
A A A和 H H H求解出来以后,就可以根据 ( 3 ) (3) (3)利用矩阵转换,推导得到外参 R R R和 t t t。
5. 利用解算出的内外参作为初值,考虑畸变参数进一步优化
以上的计算步骤并没有考虑畸变的问题,也即是默认没有畸变。但是因为有畸变,所以以上解算的结果只能作为近似结果,作为进一步优化的初值。
畸变是这么考虑的:
{
x
d
i
s
t
=
x
+
x
(
k
1
r
2
+
k
2
r
4
+
k
3
r
6
)
+
2
p
1
x
y
+
p
2
(
r
2
+
2
x
2
)
y
d
i
s
t
=
y
+
y
(
k
1
r
2
+
k
2
r
4
+
k
3
r
6
)
+
2
p
2
x
y
+
p
1
(
r
2
+
2
y
2
)
\left\{\begin{matrix} x_{dist}=x+x(k_1r^2+k_2r^4+k_3r^6)+2p_1xy+p_2(r^2+2x^2)\\ y_{dist}=y+y(k_1r^2+k_2r^4+k_3r^6)+2p_2xy+p_1(r^2+2y^2) \end{matrix}\right.
{xdist=x+x(k1r2+k2r4+k3r6)+2p1xy+p2(r2+2x2)ydist=y+y(k1r2+k2r4+k3r6)+2p2xy+p1(r2+2y2)
其中,
k
1
,
k
2
,
k
3
k_1,k_2,k_3
k1,k2,k3和
p
1
,
p
2
p_1,p_2
p1,p2分别表示径向和切向畸变;
x
,
y
x,y
x,y表示对应的归一化后的平面坐标
(
疑
问
1
)
\color{#FF0000}{(疑问1)}
(疑问1)。
x
d
i
s
t
,
y
d
i
s
t
x_{dist}, y_{dist}
xdist,ydist是像空间坐标系下面的坐标,也即是世界坐标系经过外参转换之后的坐标。畸变模型即是在(归一化后的)像空间坐标系加入了一个新的、由畸变参数描述的变换。
进一步的,经由内参矩阵得到:
{
μ
=
f
x
x
d
i
s
t
+
μ
0
ν
=
f
y
y
d
i
s
t
+
ν
0
(5)
\left\{\begin{matrix} \mu=f_xx_{dist}+\mu_0\\ \nu=f_yy_{dist}+\nu_0 \end{matrix}\right. \tag{5}
{μ=fxxdist+μ0ν=fyydist+ν0(5)
这样就可以构造重投影误差了。
构建最优化目标函数,最小化重投影误差:
∑
i
=
1
n
∑
j
=
1
m
∥
m
^
i
j
−
m
^
i
j
′
(
A
,
k
1
,
k
2
,
k
3
,
p
1
,
p
2
,
R
i
,
t
i
,
M
i
j
)
)
∥
2
\sum_{i=1}^{n}\sum_{j=1}^{m}\left \|\hat m _{ij}-\hat m_{ij} ^{'}(A, k_1, k_2, k_3, p_1, p_2, R_i, t_i, M_{ij})) \right \|^2
i=1∑nj=1∑m∥∥∥m^ij−m^ij′(A,k1,k2,k3,p1,p2,Ri,ti,Mij))∥∥∥2
其中,
m
,
n
m,n
m,n分别代表角点的个数和图片的个数。注意每张图像的外参都是不同的,但是相机的内参和畸变参数是固定的。
这个公式写的可能不太规范,减号的坐标表示的是实际的量测像点坐标(像素为单位),减号的后面是物点坐标转换来的,即通过式
(
5
)
(5)
(5)。
迭代优化后,最后的解算出来的参数包括相机的内参和畸变参数,以及每一张图像的外参。
补充知识
f x , f y f_x,f_y fx,fy的单位是像素
f f f的单位式米, f x = f d x f_x=\frac{f}{d_x} fx=dxf,其中 d x d_x dx表示x方向每个像素占多大,即一个像素在相机感光板上的物理长度。
基本的公式推导以及与共线条件方程的对应关系
相机坐标系->像平面坐标->像素坐标:
[
μ
ν
1
]
=
1
z
c
[
f
x
γ
μ
0
0
f
y
ν
0
0
0
1
]
[
x
c
y
c
z
c
]
(6)
\begin{bmatrix} \mu\\ \nu\\ 1 \end{bmatrix}= \frac{1}{z_c}\begin{bmatrix} f_x& \gamma & \mu_0\\ 0& f_y & \nu_0 \\ 0& 0& 1 \end{bmatrix}\begin{bmatrix} x_c\\ y_c\\ z_c \end{bmatrix} \tag{6}
⎣⎡μν1⎦⎤=zc1⎣⎡fx00γfy0μ0ν01⎦⎤⎣⎡xcyczc⎦⎤(6)
这个回答了疑问1,所谓归一化的坐标,就是在像空间(相机坐标系)把平面两个方向都除以Z方向。 ( 1 ) (1) (1)和 ( 2 ) (2) (2)中 z c z_c zc就是这么来的。
这里的“归一化”,对应着《摄影测量》里面共线条件方程推导中的“共线”:
考虑地更全面一点,相机标定一共涉及四个坐标系:世界坐标系、相机坐标系、图像坐标系、像素坐标系,总的转换关系为:
z
c
[
μ
ν
1
]
=
[
1
d
x
−
c
o
t
(
θ
)
d
x
μ
0
0
1
d
y
s
i
n
(
θ
)
ν
0
0
0
1
]
[
f
0
0
0
0
f
0
0
0
0
1
0
]
[
R
T
0
1
]
[
X
Y
Z
1
]
z_c\begin{bmatrix} \mu\\ \nu\\ 1 \end{bmatrix}=\begin{bmatrix} \frac{1}{d_x}& \frac{-cot(\theta)}{d_x} & \mu_0\\ 0& \frac{1}{d_ysin(\theta)} & \nu_0 \\ 0& 0& 1 \end{bmatrix}\begin{bmatrix} f&0 &0&0\\ 0&f &0&0 \\ 0& 0& 1& 0 \end{bmatrix}\begin{bmatrix} R&T\\ 0&1 \end{bmatrix}\begin{bmatrix} X\\ Y\\ Z\\ 1 \end{bmatrix}
zc⎣⎡μν1⎦⎤=⎣⎢⎡dx100dx−cot(θ)dysin(θ)10μ0ν01⎦⎥⎤⎣⎡f000f0001000⎦⎤[R0T1]⎣⎢⎢⎡XYZ1⎦⎥⎥⎤
以往看公式的时候总是会困惑上式中的
z
c
z_c
zc是怎么来的:他是
[
R
T
0
1
]
[
X
Y
Z
1
]
\begin{bmatrix} R&T\\ 0&1 \end{bmatrix}\begin{bmatrix} X\\ Y\\ Z\\ 1 \end{bmatrix}
[R0T1]⎣⎢⎢⎡XYZ1⎦⎥⎥⎤变换以后的第三个轴的坐标。
这个推导的方式不如“共线条件方程”直观,但是结果是一样的。“共线条件方程”的形式为:
这对应着
(
6
)
(6)
(6)转换后的形式:
{
μ
=
f
x
x
c
z
c
+
μ
0
ν
=
f
y
y
c
z
c
+
ν
0
\left\{\begin{matrix} \mu=f_x\dfrac{x_{c}}{z_c}+\mu_0\\ \nu=f_y\dfrac{y_{c}}{z_c}+\nu_0 \end{matrix}\right.
⎩⎨⎧μ=fxzcxc+μ0ν=fyzcyc+ν0
不管哪种方式的推导,注意其中的像方的坐标以及 f f f的单位,以及坐标的正负方向要一致。
疑问
-
拍摄过程中,标定板能不能移动位置?
拍摄标定板的时候,我以前理解的是保持标定板不动,改变相机的位置和姿态进行标定,因为这个时候物方是固定的,我只要有每张照片的角点们的坐标就可以按照张正友方法,计算单行性矩阵,然后B然后K,最后解算得到了内参,顺便也得到了每张照片的外参。得到了初值后,也可以考虑畸变通过最小化重投影误差进行整体优化。
但是似乎,网上的资料,并不强调标定板不动,我就不理解了,如果标定板也在动,相机也在动,那物方的控制怎么得到呢。回答:因为单应性矩阵是能够也确实是各张照片单独求解的,所以后面联合多片求解的时候,求解得是各片共享的内参,而且这个过程中,用到的是H,而不是什么物方的控制。所以各次拍摄中,标定板移不移动都可以解算。我会出现这个疑问是因为,思路被摄影测量中学习到的后方交会给限制了,总觉得要有物方的控制点(固定不变的已知坐标的地面点)才能够解算出未知数。
-
怎么验证校正的精度?
一般重投影误差很小的话,标定结果均可用。但很多情况下,重投影误差小,并不能说明标定参数好,要根据应用来评价的,比如目标定位精度。 -
“至于什么样的标定图像,即棋盘格怎么摆放,角点密度怎么样,棋盘格在画幅中所占比例等等,能够得到比较好的标定参数,也就是什么样的图片时有效图片,楼主可以再去试试,目前没有结论。” 真的是这样嘛?实践出真知?
-
标 定 结 果 跟 景 深 的 关 系 是 什 么 呢 ? 地 面 标 定 后 用 于 航 拍 , 效 果 如 何 ? \color{#FF0000}{标定结果跟景深的关系是什么呢?地面标定后用于航拍,效果如何?} 标定结果跟景深的关系是什么呢?地面标定后用于航拍,效果如何?
这个问题的分析见下一篇博客。
代码
这里有matlab版本的标定方法,似乎按照这个方法就可以进行啦。还有这里:
https://blog.csdn.net/JennyBi/article/details/85764988
python版本:
https://github.com/1368069096/Calibration_ZhangZhengyou_Method
C++版本:
https://blog.csdn.net/zhaitianyong/article/details/111059716