张正友标定算法学习笔记-0

张正友标定算法学习笔记

常见的相机标定方法包括基于控制场的标定方法、基于自检校模型的光束法平差、及基于二维平面格网的标定方法等。

自标定方法在不需建立标定场的情况下,利用目标周围的影像与影像之间的对应关系对相机进行标定,该方法较为灵活且成本低,缺点是稳健性较差。比较有代表性的是张正友标定法。

基本流程

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=ATA1矩阵,并利用 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=1nj=1mm^ijm^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=zc1fx00γfy0μ0ν01xcyczc(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=dx100dxcot(θ)dysin(θ)10μ0ν01f000f0001000[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

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值