在机器视觉领域，相机标定是不可或缺的一步，而张氏标定法，因其灵活性，是各种标定方法中最有影响力的一种，称其为相机标定的 "奥林匹斯山" 并不为过
张正友 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 目标函数
参考
卓越成像的11条最佳实践原则
Accurate Camera Calibration using Iterative Refinement of Control Points，2009
Calibration Best Practices
计算机与机器视觉中的高精度相机标定 (下)
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
