标定学习笔记（十）-- 一种使用校正模板的非线性相机标定方法

本文链接：https://blog.youkuaiyun.com/qq_45006390/article/details/122194824

本文归纳了郑俊，龚声蓉，刘纯平学者，于2011年发表的《一种使用矫正模板的非线性摄像机标定方法》，重在厘清相关思路，仅供学习。

1 自标定与标定模板标定

通过使用图像特征，如直线等几何不变性，并基于场景中的直线经过透视投影变换后所成的图像仍是直线的事实实现畸变校正；与通过不同视角拍摄的多幅图像之间的点的对应关系来求取镜头畸变参数。由于非可量测标定或自标定方法从本质上来讲都是基于绝对二次曲线或绝对二次曲面的方法，其最大的不足是其鲁棒性差，一般都要通过求解复杂的非线性方程来计算摄像机的内参数，从而导致运算速度和结果的精度都不理想。

可量测的非线性摄像机镜头畸变校正方法（下简称为RMC方法），可以被看成是分为两个过程：畸变模板校正和摄像机标定，然后被校正的模板图像将被用于后一步的摄像机标定中。该方法选择使用共线的 4 个点作为基准点，首先利用图像上的点满足的直线方程以及共线点的交比性质对畸变点进行校正，从而达到对畸变模板进行校正的目的，然后再完成摄像机的标定。该方法有较高的校正精度，但是，实际上只有无畸变的情况下才能保证选择的这些点是完全共线的，因此，需要对这些可能发生畸变的点进行校正。

2 RMC畸变校正方法分析

RMC方法虽然消除了非可量测标定或自标定方法鲁棒性差、运算结果精度不高等问题，但同时也导致了另外一些问题的存在，这里主要是算法复杂度的提升。该方法需要计算共线点的交比值，而由交比的性质可知，每计算一次交比就要求取若干距离和距离的比值，因此计算量将会很大，会直接影响算法的效率。

3 改进的方法

采用 2 个相邻点作为基准点，利用图像上的点满足的直线方程以及相邻点之间的距离相等的性质来完成对模板图像的校正，进而完成摄像机标定。

3.1 模板图像校正

3.1.1 校正图像点使其满足特定的直线方程

摄像机镜头畸变的存在，导致标定模板图像上的一些点会发生一定的位置变动，因此便不再位于原来的直线上。而标准模板图像（未发生畸变）上的点都在其应该在的位置上，为例提高后续实验的精度需要将这些发生畸变的点校正到其原来的位置上去。

做出如下约定：将模板图像上的 $m \times n$ 个点按照从左至右、从上往下的顺序进行编号，分别表示为 $q_1,q_2,...,q_m,q_{m+1},...,q_{2m},q_{2m+1},...,q_{n \times m}$ 。

如果点 $q_{i} = (x_i,y_i)$ 在直线 $l$ 上，则必有如下直线方程(1)成立：

$a_lx_i + b_ly_i + c_l = 0$

在这里， $a_l,b_l,c_l$ 表示定义直线 $l$ 的参数，这些参数可以利用直线上的 $m$ 个点（其中一组为 $q_1,q_2,...,q_m$ ）通过最小化方程式(2)来实现：

$P_1 = \sum_{l=1}^{n}\sum_{i =1}^{m}(a_lx_{i+(l-1)\cdot m} + b_ly_{i+(l-1)\cdot m} + c_l)^2$

由于所选用的模板有位于正交方向上的两组直线，且每一个点都是其中两条直线的交点，当上式取得最小化时，校正后的点会同时满足两条直线的方程。

3.1.2 校正图像点使相邻点之间的距离相等

在这一部分的校正过程中RMC方法利用的是交比不变性的原理，即当模板图像不存在畸变时任意共线的 4 个点所成的交比值都是相等的，且都等于 $CR(A',B',C',D')$ ，而当图像发生了畸变后，该方法通过将下面的方程式(3)最小化来实现点的畸变校正：

$J_{CR}=\sum_{l=1}^{n}\sum_{k=1}^{m-3}(CR(q_{k,l},q_{k+1,l},q_{k+2,l},q_{k+3,l})-CR(A',B',C',D'))^2$

在这里， $n$ 是模板图像上的直线个数， $m$ 是每条直线上的点数， $q_{k,l}$ 是直线 $l$ 上的点 $k$ 。如果使用这样的方法那么计算量将会非常大，为此在这里利用正方形边长相等（即相邻点之间的距离都相等）的原理来实现畸变点的校正。

定理在畸变点被校正为满足特定直线方程的前提下，若图像上所有相邻点之间的距离相等，且都等于同一个值，则模型图像将得到校正。

由于所使用的标定模板是由诸多小正方形单元方格组成的，棋盘格上相邻的两个点（像 $A,B$ 这样的两个点）之间的距离都应当是相等的，且都等于一个值（小正方形单元放个的边长）。当有畸变存在时，相邻两个点之间的距离就可能会发生变化，为例提高后续实验的精度，需要将畸变点校正到其原来的位置上去，也即是对模板图像进行校正。

不失一般性，可设直线 $l$ 上有一点 $q_i=(x_i,y_i)$ ，则在相互正交的方向上其相邻点可分别设为 $q_{i+1}=(x_{i+1},y_{i+1})$ 和 $q_{i+m}=(x_{i+m},y_{i+m})$ ，当不存在畸变时，根据正方形单元方格边长相等的，必有方程(4)(5)成立：

$\sqrt{(x_i-x_{i+1})^2+(y_i-y_{i+1})^2} = L$

$\sqrt{(x_i-x_{i+m})^2+(y_i-y_{i+m})^2} = L$

在这里， $L$ 为小正方形单元格的边长，它可以由标准模板得到。

当有畸变存在时，图像上相邻两点之间的关系将不能同上满足上面两个方程，这时需要对畸变点进行校正，可以通过最小化方程式(6)来实现：

$P_2 = |\sum_{l=1}^{n-1}\sum_{i=1}^{m-1}(\sqrt{(x_{i+(l-1)\cdot m}-x_{i+(l-1)\cdot m+1})^2+(y_{i+(l-1)\cdot m}-y_{i + (l-1)\cdot m +1})^2 } + \sqrt{(x_{i+(l-1)\cdot m} - x_{i+l \cdot m})^2+(y_{i+(l-1)\cdot m}-y_{i+l \cdot m})^2} -2L|$

考虑到边界处的特殊性，将上式转换为方程式(7)：

$P_2' = |\sum_{l=1}^{n-1}\sum_{i=1}^{m-1}({(x_{i+(l-1)\cdot m}-x_{i+(l-1)\cdot m+1})^2+(y_{i+(l-1)\cdot m}-y_{i + (l-1)\cdot m +1})^2 } + {(x_{i+(l-1)\cdot m} - x_{i+l \cdot m})^2+(y_{i+(l-1)\cdot m}-y_{i+l \cdot m})^2} -2L^2|$

考虑到边界处的特殊性，在 $i=m$ 和 $l = n$ 情况下，可分别得到下面两个方程(8)(9)：

$P_3 = |\sum_{l=1}^{n-1}((x_{l \cdot m} - x_{(l+1)\cdot m})^2 + (y_{l \cdot m}-y_{(l+1)\cdot m} )^2 - L^2)|$

$P_4 = |\sum_{i=1}^{m-1}((x_{i +(n-1)\cdot m} - x_{i + (n-1)\cdot m+1})^2 + (y_{i +(n-1)\cdot m}-y_{i+(n-1)\cdot m+1} )^2 - L^2)|$

3.1.3 标定模板图像的校正

上面两小节对标定用的模板图像进行分析，当不存在畸变时模板图像上的点应满足的条件（式(1)、(4)、(5)所示）以及存在畸变时将畸变点校正为理想点时应满足的条件(式(2)、(7)、(8)、(9)所示)。为了将已经发生畸变的模板图像校正为标准模板图像，需要将畸变点校正为理想点（畸变点本来所在的位置），则需要将畸变点校正为理想点满足的所有相关约束联立起来，将整体的误差限定到最小。

综合式(2)、(7)、(8)、(9)，得到下式(10)：

$P = P_1 + P_2' +P_3 + P_4$

至此畸变点校正的问题就转为求取 $min(P)$ （函数 $P$ 的最小值），当图像点都位于相应的直线上且相邻两点之间的距离都等于一固定值（小正方形单元方格的边长）时，式(10)将会取得最小值。

这里给出相关描述：通过使用迭代非线性的 Levenberg-Marquardt(LM)算法来实现最小化过程，以小正方形单元格的边长 $L$ 、畸变点 $q_{i,d}$ 的集合和初始估计的直线参数 $(a_l,b_l,c_l)$ 作为输入参数。

$L$ 是一个固定值，可以从模板直接得到；畸变点 $q_{i,d}$ 可以从图像中提取出来；直线参数 $(a_l,b_l,c_l)$ 的初始估计值可以由畸变的图像点 $q_{i,d}$ 得到。在每一次的迭代过程中，畸变点 $q_{i,d}$ 的位置和直线参数 $(a_l,b_l,c_l)$ 都会被重新计算一次；当畸变点 $q_{i,d}$ 的位置逐渐接近于理想点位置 $q_{i,p}$ ，此时 $q_{i,p}$ 满足由参数 $(a_l,b_l,c_l)$ 确定的直线且相邻点之间的距离无限接近于 $L$ ，此时非线性搜索将会结束，式(10)即可取得最小值。

3.2 非线性摄像机镜头畸变模型建立

由于摄像机制造工艺上的不足，摄像机镜头的成像过程往往不服从针孔模型，也即非线性模型。常见的非线性镜头畸变有三种类型，分别是径向畸变、切向畸变和薄棱镜畸变。

图像的径向畸变是指图像点从它的理想位置向内或向外移动，即在径向上出现误差；而图像的切向畸变是指由于装配误差，组成光学系统的光学镜头的光轴不可能完全共线，即图像点在切向上出现了偏差；薄棱镜畸变是指由光学镜头制造误差和成像敏感阵列制造误差引起的图像变形。

常见的非线性镜头畸变模型是包含上述三种畸变类型的多项式畸变模型（Polynomial distortion model，NLPD），该畸变模型定义了一个从畸变点 $q_d = (u_d,v_d)$ 到理想点 $q_p = (u_p,v_p)$ 之间的映射，该映射可以描述如下四个方程(11)(12)(13)(14)：

$u_d = u_p + \delta_u$

$v_d = v_p + \delta_v$

$\delta_u = \Delta u_d \cdot (k_1 \cdot r_d^2 + k_2 \cdot r_d^4 + ...) + p_1 \cdot (3 \Delta u_d^2 + \Delta v_d^2)+2 p_2 \cdot \Delta u_d \cdot \Delta v_d + s_1 \cdot r_d^2$

$\delta_v = \Delta v_d \cdot (k_1 \cdot r_d^2 + k_2 \cdot r_d^4 + ...) + p_2 \cdot (3 \Delta v_d^2 + \Delta u_d^2)+2 p_1 \cdot \Delta u_d \cdot \Delta v_d + s_2 \cdot r_d^2$

这里定义图像的畸变中心为 $c=(u_0,v_0)$ ，则有 $\Delta u_d = u_d - u_0$ ， $\Delta v_d = v_d - v_0$ ， $r_d^2 = \Delta u_d^2 + \Delta v_d^2$ 。 $k_1,k_2$ 等是径向畸变系数， $p_1,p_2$ 是切向畸变系数， $s_1,s_2$ 是薄棱镜畸变系数。

针对三种畸变类型之间的关系，其中径向畸变（主要是前两阶的）和切向畸变起主导作用，而薄棱镜畸变只是会引起额外的径向畸变和切向畸变，而且考虑过多的畸变会导致结果的不稳定，薄棱镜畸变往往忽略不计。因此，可以将NLPD的式(13)和(14)简化为下面两式(15)(16)：

$\delta_u = \Delta u_d \cdot (k_1 \cdot r_d^2 + k_2 \cdot r_d^4) + p_1 \cdot (3 \Delta u_d^2 + \Delta v_d^2) + 2p_2 \cdot \Delta u_d \cdot \Delta v_d$

$\delta v = \Delta v_d \cdot (k_1 \cdot r_d^2 + k_2 \cdot r_d^4) + p_2 \cdot (\Delta u_d^2 + 3 \Delta v_d^2) + 2 p_1 \cdot \Delta u_d \cdot \Delta v_d$

3.3 摄像机标定

一般来说，摄像机标定可以分为两个大类：传统的摄像机标定方法呵摄像机自标定方法。传统的摄像机标定方法标定精度比较高，但是标定过程复杂，需要使用结构已知、精度很高的标定块作为空间参照物；自标定方法不需要使用标定块，但精度相对较低，鲁棒性比较差。张正友标定法实际上是传统标定方法和自标定方法的一个折中，克服了两种方法的确定，拥有较高的标定精度，且鲁棒性较好。

张的方法仅考虑了径向畸变，忽略了切向畸变，在这里，采用式(15)(16)所确定的包含径向畸变和切向畸变的两种类型畸变的摄像机镜头畸变模型进行摄像机标定。根据张的方法可以求出除切向畸变系数 $p_1,p_2$ 之外的所有参数，这里给出 $p_1,p_2$ 的求解。

式(15)和式(16)给出了图像上的理想点和畸变点之间的一种映射关系，根据张的方法求径向畸变系数 $k_1,k_2$ 的原理，这里将 $k_1,k_2,p_1,p_2$ 放在一起进行求解，下面给出式(15)和(16)的一个变形式(17)：

$\begin{bmatrix} \Delta u_d \cdot (\Delta u_d^2+ \Delta v_d^2) & \Delta u_d \cdot (\Delta u_d^2+ \Delta v_d^2)^2 & 3\Delta u_d^2+\Delta v_d^2 & 2\Delta u_d \cdot \Delta v_d \\ \Delta v_d \cdot (\Delta u_d^2+ \Delta v_d^2) & \Delta v_d \cdot (\Delta u_d^2+ \Delta v_d^2) & 2 \Delta u_d \cdot \Delta v_d & \Delta u_d^2 + 3\Delta v_d^2 \end{} \cdot \begin{bmatrix} k_1 \\ k_2 \\ p_1 \\ p_2 \end{} = \begin{bmatrix} \delta u \\ \delta v \end{}$