张正友标定法

大二下学期只是听明白了方法,却搞不懂究竟怎么算的

单应矩阵

单应矩阵给出从一个平面到另一个平面的映射,它包括摄像机内外参数矩阵

张正友标定概览

前提

  • 内参数矩阵:五参数模型
  • 畸变模型:四阶径向畸变模型
  • 标定物:平面靶标
  • 世界坐标系W按照下图置于靶标平面,原点设在靶标一角,$X_w$和$Y_w$方向沿着靶标平面,$Z_w$方向垂直于靶标平面($z=0$平面上)(2011考点,设置方法和理由)

基本步骤

  1. 不考虑畸变,标定摄像机参数,得到参数的线性初值
  2. 利用线性初值,进行非线性标定。得到畸变参数
  3. 重复1和2,直至参数收敛

具体步骤

  1. 求解单应矩阵
  2. 求解摄像机内参数
  3. 求解摄像机外参数
  4. 求解畸变系数$k_1$, $k_2$
  5. 参数最优化

求解单应矩阵

单应矩阵的推导

显而易见,所求的单应矩阵是从世界坐标系的棋盘格平面到摄像机图像平面的映射,而这与摄像机内外参数有着直接的联系。
设靶标上一点在世界坐标系的位置为$[x_{wi} \quad y_{wi} \quad z_{wi}]$,其在图像上对应的坐标是$[u_i \quad v_i]$,图像与空间点的映射关系为

由于$z_{wi}=0$,可将上式化简为

其中3x3矩阵 $H=M_{in}[r_1 \quad r_2 \quad P]$ 是单应矩阵,定义$H=[h_1 \quad h_2 \quad h_3]$ (2011考点,单应矩阵的形式和物理意义)
(为了码个矩阵还费了好大劲,改了mathjax的解析包,ps. npm命令突然不好使卸载重装node.js就ok了)
由于$r_1$和$r_2$是旋转矩阵$R$的单位正交列向量,有$r_1^T r_2=0$,$|r_1|=|r_2|=1$,得到以下两个约束方程

求解步骤

论文中给出了一种最大似然估计的方法,考虑零均值高斯噪声,协方差矩阵为$\Lambda_i=\sigma^2 \textbf{I}$
设$H=\begin{bmatrix} \bar{h_1^T}\\ \bar{h_2^T}\\ \bar{h_3^T}\\ \end{bmatrix}$,9维列向量$x=\begin{bmatrix}\bar{h_1} \\ \bar{h_2}\\ \bar{h_3}\\ \end{bmatrix}$,设$F_i$和$f_i$分别为物体坐标和图像坐标
对于每一对对应点,成像方程为下式

若有N对对应点,将有$Lx=0$,其中L是2Nx9的矩阵,通过非线性最小二乘法求出x,进而得到H
H共八个未知量,至少提取4个角点可解

求解内参数矩阵

设$B$为以下3x3对称矩阵

其中$M_{in}=\begin{bmatrix} k_x&k_s&u\\0&k_y&v\\0&0&1\\ \end{bmatrix}$
$B_{11}=\frac{1}{k_x^2}$
$B_{12}=B_{21}=-\frac{k_s}{k_x^2 k_y}$
$B_{13}=B_{31}=\frac{k_s v-k_y u}{k_x^2 k_y}$
$B_{23}=B_{32}=\frac{k_s(k_y u-k_s v)}{k_x^2 k_y^2}-\frac{v}{k_y^2}$
$B_{33}=\frac{(k_s v-k_y u)^2}{k_x^2 k_y^2}+\frac{v^2}{k_y^2}+1$

因此可以由一个六维列向量唯一确定$\textbf{b}=\begin{bmatrix}B_{11}&B_{12}&B_{22}&B_{13}&B_{23}&B_{23} \end{bmatrix}^T$
可推得

令$v_{ij}=\begin{bmatrix} h_{1i}h_{1j}&h_{1i}h_{2j}+h_{2i}h_{1j}&h_{2i}h_{2j}&h_{1i}h_{3j}+h_{3i}h_{1j}&h_{2i}h_{3j}+h_{3i}h_{2j}&h_{3i}h_{3j} \end{bmatrix}^T$
可得到$h_i^T B h_j=v_{ij}^T \textbf{b}$
由3.1.1的两个约束方程可得到如下方程

获取n张不同的靶标图像,得到方程$\textbf{Vb}=\textbf{0}$,其中$\textbf{V}$是2nx6维度的矩阵,当$n \geq3$时,六个方程可唯一解出$\textbf{b}$的基础解系,因此确定了B矩阵(带倍乘系数)至少获取三张图像
然后就能通过B矩阵计算摄像机内参数
$v=(B_{12}B_{13}-B_{11}B_{23})(B_{11}B_{22}-B_{12}^2)$
$\lambda=B_{33}-[B_{13}^2+v(B_{12}B_{13}-B_{11}B_{23})]/B_{11}$ \quad (系数)
$k_x=\sqrt{\frac{c}{B_{11}}}$
$k_y=\sqrt{\frac{\lambda B_{11}}{B_{11} B_{12}-B_{12}^2}}$
$k_s=-\frac{B_{12} k_x^2 k_y}{\lambda}$
$u=\frac{k_s v}{k_y}-\frac{B_{13} k_x^2}{\lambda}$

求解外参数矩阵

既然已经求出了内参数矩阵$M_{in}$,就可以通过单应矩阵的定义反解出各个外参数
$r_1=\alpha M_{in}^{-1} h_1, \quad r_2=\alpha M_{in}^{-1} h_2, \quad r_3=r_1 \times r_2, \quad P=\alpha M_{in}^{-1} h_3$
其中$\alpha=1/|M_{in}^{-1} h_1|=1/|M_{in}^{-1} h_2|$(ps.这个常数因子哪里来的我很费解)
由于数据必然含有噪声,因此直接计算出的$R=[r_1 \quad r_2 \quad r_3]$一般不满足旋转矩阵的性质条件,最优的旋转矩阵可以通过奇异值分解(SVD)的方法获取

求解畸变系数

论文里貌似没有这一步,可能在其他论文里有吧。得到了内外参数后,根据假设的畸变模型计算畸变系数$Ak=Q$,只需解$(x_i,y_i)$代表角点在世界坐标系下的坐标,$(u_0,v_0)$代表角点在实际图像中的坐标,$(u_1,v_1)$代表角点通过计算得到的在图像中的理想坐标。

然后通过最小二乘法计算$k=(A^TA)^{-1}A^TQ$

参数最优化

根据目标函数F,通过最大似然估计法实现参数最优化(或者叫做非线性最小二乘法也可以?),只需最小化F

其中$p_{ij}$表示世界坐标点,$I_{ij}$表示第j幅图像中第i个角点的实际坐标,$\bar{I_{ij}}$表示根据参数计算出的图像坐标,论文中没有k1和k2作为参数,不过实质是一样的,利用m幅靶标图像中的共mxn个角点来最小化目标函数

OpenCV函数

findChessboardCorners:检测标定板角点
cornerSubPix:亚像素精度角点检测
drawChessboardCorner:绘制角点

参考文献:
Zhengyou Zhang, ``A Flexible New Technique for Camera Calibration’’, IEEE Transaction on Pattern Ananlysis and Machine Intelligence, 2000.

Thanks!