翻越相机标定的奥林匹斯
在机器视觉领域，相机标定是不可或缺的一步，而张氏标定法，因其灵活性，是各种标定方法中最有影响力的一种，称其为相机标定的 "奥林匹斯山" 并不为过
张正友 99年发表的论文 "Flexible Camera Calibration by Viewing a Plane from Unknown Orientations"，在 2013年获得 ICCV 的 Helmholtz Prize，便是对其影响的认可
鉴于该论文的广泛影响，张后来再三完善论文细节，治学严谨可见一斑。网上下载最多的是 "A Flexible New Technique for Camera Calibration"，和获奖那篇差异不大
关于张氏标定法，不再赘述，本文将围绕"标定精度"，和读者一起翻越相机标定的"奥林匹斯山"
1 相机标定
1.1 定义
已知：世界坐标系中(建在平面标定板上)，几组特征点的空间坐标，以及在像素坐标系中，特征点对应的像素位置坐标
求解：相机的内参和畸变系数
标定板上特征点的空间坐标(3d)，通过相机模型，与特征点的像素位置坐标(2d)关联起来，如下：
$\qquad s \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} r_{11} & r_{12} & r_{13} & t_1 \\ r_{21} & r_{22} & r_{23} & t_2 \\ r_{31} & r_{32} & r_{33} & t_3 \end{bmatrix} \begin{bmatrix} X_w \\ Y_w \\ Z_w\\ 1 \end{bmatrix} $
1.2 二维展开
相机标定的过程，就是输入几组已知的 3d 特征点坐标和提取的 2d 特征点坐标，将反投影误差构建为目标函数，在已知相机模型方程的基础上，通过最小化目标函数，得到相机模型方程的参数
开个脑洞，类比《三体》中的质子展开过程：高维度的相机标定问题，如果展开到二维空间，就是一个曲线拟合的问题，如下图，参见 Ceres Solver 的曲线拟合例程
分析相机标定的过程，得出影响标定(精度)的三个因素：1) 特征提取；2) 最优化方法；3) 相机模型
2 特征提取
以视场(对角线) 100mm 的机器视觉系统为例，根据卓越成像的最佳实践原则 #3，可知镜头的工作距离约为 200~400mm
选 1/3 英寸(对角线6mm) CIS芯片，同样根据最佳实践原则 #3，推测出镜头的理想焦距为 12~24mm
假定 CIS芯片的分辨率为 1024(H) x 768(V)，则 H 方向的"像素分辨率"为 80mm/1024pixel ≈ 0.08mm/pixel
这意味着，如果提取的特征点偏差 1个像素，造成的尺寸偏差为 0.08mm，因此，特征提取的精度非常重要
2.1 标定板图案
特征提取，属于图像处理范畴，处理的是标定板在相机中的成像图片，提取的是标定板图案的特征点 (如角点、圆心等)
常用标定板的图案，有棋盘格、圆、非对称圆、圆环，以及 ChArUco 等，Halcon 公司有特定图案的 Halcon 标定板
一般而言，圆环标定板的标定精度最高，圆次之，棋盘格最低
2.1.1 特征提取精度
参考[2]中，对三种图案标定板的各五张合成图像(非相机拍摄)，添加相同的径向畸变，比较不同图案的特征提取精度
1) 棋盘格，特征为角点，利用 OpenCV 先 findChessboardCorners() 粗定位，再 cornerSubPix() 精定位，得定位误差 6pixel
2) 圆标定板，特征为圆心，定位误差次小，为 2.6pixel
3) 圆环标定板，特征为圆环中心，定位误差最小，为 1.7pixel
对于0.08mm/pixel 的机器视觉系统，当镜头确定时(畸变系数固定)，不同标定板对应的特征，提取精度最大相差 4.3pixel，约 0.34mm
2.1.2 反投影误差
用这三种图案的各五张合成图像，来进行标定，得到反投影误差的 RMSE 分别为 0.1263，0.0517 和 0.0515
而用相机拍摄的三种图案的实际图像，进行标定后，得到反投影误差的 RMSE 分别为 0.139，0.135 和 0.115
综上，圆环标定板的标定精度，略高于圆标定板，二者均高于棋盘格标定板
2.2 迭代标定法
特征提取出现偏差，是因为在标定过程中，标定板要变换位姿，不同位姿会引起透视收缩，导致特征检测精度降低，张正友的论文里也提到过
2.2.1 偏差原因
当棋盘格旋转一定角度时，cornerSubPix() 精定位算法中的梯度，并不和边缘垂直，导致提取的角点有偏差
将标定板进行透视校正，使棋盘格正对相机，此时算法中的梯度和图像边缘垂直，角点提取无偏差
针对位姿导致的特征提取偏差，参考[2]提出了一种迭代标定法：先将图像转换为平行正对，再检测特征，重新标定，迭代直至收敛
2.2.2 方法步骤
1）给定图像的特征检测(角点、圆心或圆环中心)
2）相机标定，得到标定参数
重复以下步骤，直至收敛
3）畸变校正和透视校正(转换为平行正对图像)
4）在平行正对图像中，检测特征(模板匹配算法)
5）根据标定参数(每次迭代更新)，将特征点转换回原始位置
6）相机标定，得到新的标定参数
选相机真实拍摄的图像，采用该方法标定，得到反投影误差的 RMSE，如下：
从表中看，采用迭代标定法，棋盘格的 RMSE 从 0.14 降到了 0.08，而圆和圆环的，则分别从 0.14 和 0.12 降到了 0.07 和 0.06
以上结果，是参考[2]中所列，本人未实践过，但在一家公司公众号的文章中 (参考[4])，发现了类似迭代标定法的动图，如下：
3 优化方法
张的论文中，输入图像 ≥11幅时，标定误差显著减小，而 Bouguet 标定包，用 20~25幅图像，图像数量越多，意味着提取的特征点也越多
实际中，受照明不均匀、标定板或镜头污染、提取算法等的限制，即使标定板转换为平行正对，提取的特征点也会有离群点
3.1 RANSAC 标定法
对于特征中离群点的剔除，常用的一种方法是 RANSAC：以反投影误差 $E_{reproj}$ 做阈值，小于的为内点，大于的为外点，不断迭代使选定的内点都满足 $<E_{reproj}$
3.1.1 方法步骤
1）张氏标定法，得到标定参数
2）计算所有特征点的反投影误差 $E_{reproj}$
3）选 $E_{reproj} < T_{reproj}$ 的所有内点，再次标定得到新的参数
重复步骤 2) 和 3)，直到所有的内点都满足 $E_{reproj} < T_{reproj}$ (文中 T 取 2)
4）对每一幅标定图像，计算其反投影误差 $E_{img}$，设定 RANSAC 参数初值，例如，$T_{img}=1.2 E_{img}$，最大内点数 $N=\infty$，置信度 $p=0.99$，$i=1$ 等
5）将提取的特征点，在像素位置坐标中画十字线，均分成四组 (避免四点共线)
6）四组中各选一个点，计算相机外参，综合之前的相机内参，筛选出满足 $< T_{img}$ 的内点集 $S_{in}^i$
7）如果本次筛选的 $S_{in}^i$，比之前筛选的内点数量多，则更新 $N$ 值为 $S_{in}^i$ 的内点数
8）当 $i > N$ 时，进入步骤 9)，否则，令 $i=i+1$，重复步骤 5) 6) 7）
9）对每一幅图像，都执行步骤 5) 6) 7) 8)，获取每一幅图像的最大内点集 $S_{con}$
10）综合每幅图像，选出的内点集，再次标定，得到最终的标定参数
以上步骤，和 OpenCV 中的 findHomography() 函数，在求解单应性矩阵时所用的 RANSAC 法类似，可参考之
3.1.2 测试效果
对于一幅合成的标定图像，添加高斯噪声后，采用阈值法 和 RANSAC 法，筛选出的内点集，如下：
对于一幅真实的标定图像，当因光照不均等存在离群点时，采用阈值法 和 RANSAC 法，筛选效果如下：
从统计学上分析，随着提取特征点中离群点的剔除，标定参数更逼近真实值，相机的标定精度也越来越高
由此可知，采用阈值法 和 RANSAC 法，对于大部分的离群点，都能成功的筛选出来
3.2 目标函数
3.2.1 2d 和 3d
考虑镜头的畸变，相机标定是一个非线性优化的过程，对于张氏标定法，目标函数在 2d 像平面上，如下：
可定义为像平面中检测到的特征点 $m_{u}^d$，和反投影到像平面上的特征点$\hat m_{u}^d$之间的像素差，如下：
$\quad\displaystyle\sum_{i=1}^{N} \sum_{j=1}^{L}|| m_{u,i,j}^d - \hat{m}_{u,i,j}^{d}(f_x, f_y, u_0, v_0, k_1, k_2, R_i, t_i)||^2$
在实际应用中，视觉系统测量的是 3d 空间中的特征点，2d 像平面中像素差相同的两组点对，投射到 3d 空间中的距离差反而不相同
可重定义目标函数，为检测到的特征点转换到相机坐标系 $m_c$，和反投影到相机坐标系上的特征点 $\hat m_c$ 之间的距离差，如下：
$\quad\displaystyle\sum_{i=1}^N \sum_{j=1}^L || m_{c,i,j}(R_i, t_i) - \hat{m}_{c, i, j}(f_x, f_y, u_0, v_0,k_1, k_2, R_i, t_i)||^2$
拍摄实际的图像，分别用上述两种方法进行标定，结果表明：二者的反投影误差几无差别
3.2.2 $E_1$ 和 $E_2$
目标函数变了，用原来的反投影误差来评价标定精度，不再合适，为此，引入三个新的评价指标：
$\quad E_1 = \begin{split}\frac{1}{n}\displaystyle\sum_{i=1}^{n}\sqrt{||M_{c,i} - \hat M_{c,i}||^2} \end{split}$
$\quad E_2 = \begin{split}\frac{1}{n}\displaystyle \sum_{i=1}^{n}\sqrt{||M_{c,i} - \hat L_{c,i}||^2} \end{split}$
$\quad E_3 = \begin{split}\frac{1}{m}\displaystyle \sum_{i,j=1}^{n}\sqrt{(|| M_{w,i} - M_{w,j}|| - ||\hat M_{c,i} - \hat M_{c,j}||)^2} \end{split}$
再次比较这两种标定方法，结果表明：3d 目标函数的优化精度要高于 2d 目标函数
参考
Accurate Camera Calibration using Iterative Refinement of Control Points，2009
Accurate and robust estimation of camera parameters using RANSAC，2012
Camera Calibration Toolbox for Matlab,Jean-Yves Bougue
OpenCV Tutorials / feature2d module / Basic concepts of the homography explained with code
A novel optimization method of camera parameters used for vision measurement，2013
Review of Calibration Methods for Scheimpflug Camera，2018
Automatic machine vision calibration using statistical and neural network methods，2005