目录
基本概念
通常采用如上图所示的小孔成像模型（pinhole camera model）研究摄像机标定问题。
如上图所示，摄像机标定涉及的四大坐标系：
- 世界坐标系$({X}_w, {Y}_w, {Z}_w)$：原点为$O_w$；
- 摄像机坐标系$({X}_c, {Y}_c, {Z}_c)$：原点$O_c$为光学中心，通常是小孔（pinhole）的位置；
- 像平面坐标系／图像物理坐标系$({x}, {y})$：通过光学系统在感光芯片上成像的图像坐标系；
- 像素坐标系$({u}, {v})$。
在三维图像处理中，标定就是确定摄像机的内外参数（intrinsic and/or extrinsic camera parameters）：
- 内部参数：描述摄像机的内部几何结构和光学特性；
- 外部参数：描述摄像机坐标系与世界坐标系的关系。
为了从2张或更多的图像中重构环境的三维信息，有必要知道对象的三维环境（世界坐标系）与二维图像之间的关系：
- ［三维 二维］透视投影（perspective projection）：已知摄像机的投影矩阵，可近似确定三维点的二维投影。
- ［二维 三维］逆投影（backprojection）：若已知图像上的二维点，可确定三维空间中三维点所在的射线。若存在三维点两个或更多的视图（view），其坐标可通过三角测量法（triangulation）确定。
透视投影可以减少场景分析（scene analysis）过程中处理的数据量，逆投影可用于从二维图像的特征中重建三维信息。
真实的摄像机和镜头通常都有成像误差，并不满足小孔成像模型的约束条件。成像误差主要来源于几方面：
- 成像设备分辨率低，导致空间分辨率低；
- 大多数镜头不对称，并且会导致畸变；
- 装配精度不够（比如，图像传感器的中心不在光轴上［optical axis］，传感器与镜头不平行）；
- 摄像机硬件和采集硬件之间的计时误差。
摄像机标定可在线或离线完成，可采用自标定（self-calibration）或机器学习的方法。通常的方法是借助标定物：（i）确定摄像机的参数；（ii）确定摄像机坐标系和世界坐标系之间的变换。
坐标系
一、像素坐标系与像平面坐标系
假定每个像素在$u$轴和$v$轴方向的物理尺寸分别为$\mathrm{d}x$和$\mathrm{d}y$，那么就有：
\begin{equation} \begin{aligned} u&={x\over \mathrm{d}x}+u_0\\ v&={y\over \mathrm{d}y}+v_0。 \end{aligned} \end{equation}
$\mathrm{d}x$、$\mathrm{d}y$表示每个像素在感光芯片上占据的实际大小，通过它建立了像素坐标系和真实尺寸坐标系的联系，$u_0$、$v_0$是像素平面中心。上述公式可以表示为矩阵形式
\begin{equation} \begin{bmatrix}u\\ v\\ 1\end{bmatrix} = \begin{bmatrix} {1\over \mathrm{d}x} & 0 & u_0 \\ 0 & {1\over \mathrm{d}y} & v_0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}x\\ y\\ 1\end{bmatrix}， \label{eq:pixel-vs-image-coordination} \end{equation}
或者
\begin{equation} \begin{bmatrix}x\\ y\\ 1\end{bmatrix} = \begin{bmatrix} {\mathrm{d}x} & 0 & -u_0\mathrm{d}x \\ 0 & {\mathrm{d}y} & -v_0\mathrm{d}y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}u\\ v\\ 1\end{bmatrix}。 \end{equation}
该变换建立了像素与物理尺寸的对应关系。
二、摄像机坐标系与世界坐标系
这两个坐标系的关系可用旋转矩阵$\mathbf R$和平移矩阵$\mathbf t$表示如下：
\begin{equation} \begin{bmatrix}X_c\\ Y_c\\ Z_c \\ 1\end{bmatrix} = \begin{bmatrix} \mathbf R & \mathbf t \\ \mathbf 0^\top & 1 \end{bmatrix} \begin{bmatrix}X_w\\ Y_w\\ Z_w \\ 1\end{bmatrix} = \mathbf L_w \begin{bmatrix}X_w\\ Y_w\\ Z_w \\ 1\end{bmatrix}。 \end{equation}
其中，$\mathbf R$为$3\times 3$的矩阵，$\mathbf t$为$3\times 1$的矩阵，$\mathbf 0^\top = [0\; 0\; 0]$。
三、摄像机坐标系与像平面坐标系（成像投影关系）
对线性摄像机模型（针孔模型），可得如下关系（$f$为焦距）
\begin{equation} \frac{x}{f} = \frac{X_c}{Z_c} \\ \frac{y}{f} = \frac{Y_c}{Z_c}， \end{equation}
也就是
\begin{equation} xZ_c = fX_c \\ yZ_c = fY_c， \end{equation}
同样可以表示为矩阵形式
\begin{equation} Z_c \begin{bmatrix}x\\ y\\ 1\end{bmatrix} = \begin{bmatrix} f & 0 & 0 & 0 \\ 0 & f & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix}X_c\\ Y_c\\ Z_c\\ 1\end{bmatrix}。 \end{equation}
四、像素坐标系与世界坐标系
根据上述坐标系的关系，可以推导出像素坐标系与世界坐标系的关系
\begin{equation} \begin{aligned} Z_c \begin{bmatrix}u\\ v\\ 1\end{bmatrix} &= \begin{bmatrix} {1\over \mathrm{d}x} & 0 & u_0 \\ 0 & {1\over \mathrm{d}y} & v_0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} f & 0 & 0 & 0 \\ 0 & f & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} \mathbf R & \mathbf t \\ \mathbf 0^\top & 1 \end{bmatrix} \begin{bmatrix}X_w\\ Y_w\\ Z_w \\ 1\end{bmatrix}\\ &=\begin{bmatrix} a_x & 0 & u_0 & 0 \\ 0 & a_y & v_0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} \mathbf R & \mathbf t \\ \mathbf 0^\top & 1 \end{bmatrix} \begin{bmatrix}X_w\\ Y_w\\ Z_w \\ 1\end{bmatrix} =\mathbf M_1\mathbf M_2\mathbf X_w=\mathbf M\mathbf X_w。 \end{aligned} \label{eq:pixel-vs-world-coordination} \end{equation}
$\mathbf M$是$3\times 4$的矩阵，称为投影矩阵；$\mathbf M_1$只与摄像机的内部结构有关，称为摄像机内部参数；$\mathbf M_2$（也就是$\mathbf R$和$\mathbf t$）完全由摄像机相对于世界坐标系的方位决定，称为摄像机外部参数（extrinsic parameter）。已知投影矩阵$\mathbf M$，可以根据空间点的坐标确定求出其在图像中的位置，由于$\mathbf M$不可逆，已知图像中的点，无法唯一确定其空间中的对应位置。
Faugeras标定法
直接线性变换（DLT，direct linear transformation）[1]是将像素点和物体成像几何关系在齐次坐标下写成透视投影矩阵的形式，确定该矩阵的过程称为标定。利用直接线性变换的标定法不考虑相机镜头的非线性畸变。Faugeras标定法[2, P. 56][3, P. 41]是一种直接线性变换标定方法。
公式\eqref{eq:pixel-vs-world-coordination}可以记为如下形式
\begin{equation} Z_c \begin{bmatrix}u\\ v\\ 1\end{bmatrix} = \mathbf M \begin{bmatrix}X_w\\ Y_w\\ Z_w \\ 1\end{bmatrix} = \begin{bmatrix} m_{11} & m_{12} & m_{13} & m_{14} \\ m_{21} & m_{22} & m_{23} & m_{24} \\ m_{31} & m_{32} & m_{33} & m_{34} \end{bmatrix} \begin{bmatrix}X_w\\ Y_w\\ Z_w \\ 1\end{bmatrix}， \end{equation}
展开可得
\begin{equation}
\begin{aligned}
u &= \frac{m_{11}X_w + m_{12}Y_w + m_{13}Z_w + m_{14}}{m_{31}X_w + m_{32}Y_w + m_{33}Z_w + m_{34}} \\
v &= \frac{m_{21}X_w + m_{22}Y_w + m_{23}Z_w + m_{24}}{m_{31}X_w + m_{32}Y_w + m_{33}Z_w + m_{34}}。
\end{aligned}
\end{equation}
对于无畸变（distortion-free）摄像机模型，标定就是确定参数矩阵$\mathbf M$。将上式转换为线性方程组的形式
\begin{equation} \begin{aligned} m_{11}X_w + m_{12}Y_w + m_{13}Z_w + m_{14} - m_{31}uX_w - m_{32}uY_w - m_{33}uZ_w &= m_{34}u\\ m_{21}X_w + m_{22}Y_w + m_{23}Z_w + m_{24} - m_{31}vX_w - m_{32}vY_w - m_{33}vZ_w &= m_{34}v。 \end{aligned} \end{equation}
令$m_{34}=1$（相当于方程两边都除以$m_{34}$），方程组可记为等价的矩阵形式
\begin{equation} \mathbf A\mathbf L=\mathbf U， \end{equation}
其中，
\[ \begin{aligned} \mathbf A &= \begin{bmatrix} X_{w1} & Y_{w1} & Z_{w1} & 1 & 0 & 0 & 0 & 0 & -u_1X_{w1} & -u_1Y_{w1} & -u_1Z_{w1} \\ 0 & 0 & 0 & 0 & X_{w1} & Y_{w1} & Z_{w1} & 1 & -v_1X_{w1} & -v_1Y_{w1} & -v_1Z_{w1} \\ &&&&&\ldots \\ X_{wn} & Y_{wn} & Z_{wn} & 1 & 0 & 0 & 0 & 0 & -u_nX_{wn} & -u_nY_{wn} & -u_nZ_{wn} \\ 0 & 0 & 0 & 0 & X_{wn} & Y_{wn} & Z_{wn} & 1 & -v_nX_{wn} & -v_nY_{wn} & -v_nZ_{wn} \end{bmatrix}\\ \mathbf L &= \left[m_{11}\;m_{12}\;m_{13}\;m_{14}\;m_{21}\;m_{22}\;m_{23}\;m_{24}\;m_{31}\;m_{32}\;m_{33}\right]^\top \\ \mathbf U &= \left[u_1\;v_1\ldots u_n\;v_n\right]^\top ， \end{aligned} \]
其最小二乘解为
\begin{equation} \mathbf L=\left(\mathbf A^\top\mathbf A\right)^{-1}\mathbf A^\top\mathbf U。 \label{eq:ls-solution} \end{equation}
该方程组有11个未知数，意味着使用至少6组点对（$n\geq 6$）就可以求解出全部的系数。通过$\mathbf L$得到参数矩阵$\mathbf M$后，可求解摄像机的全部内外参数。
根据\eqref{eq:pixel-vs-world-coordination}可得
\[ m_{34} \begin{bmatrix} \mathbf m_1^\top & m_{14} \\ \mathbf m_2^\top & m_{24} \\ \mathbf m_3^\top & 1 \end{bmatrix} =\begin{bmatrix} a_x & 0 & u_0 & 0 \\ 0 & a_y & v_0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} \mathbf r_1^\top & t_{x} \\ \mathbf r_2^\top & t_{y} \\ \mathbf r_3^\top & t_z\\ \mathbf 0^\top & 1 \end{bmatrix}， \]
也就是
\begin{equation} m_{34} \begin{bmatrix} \mathbf m_1^\top & m_{14} \\ \mathbf m_2^\top & m_{24} \\ \mathbf m_3^\top & 1 \end{bmatrix} =\begin{bmatrix} a_x \mathbf r_1^\top + u_0\mathbf r_3^\top & a_xt_x+u_0t_z\\ a_y \mathbf r_2^\top + v_0\mathbf r_3^\top & a_yt_y+v_0t_z\\ \mathbf r_3^\top & t_z \end{bmatrix}。 \end{equation}
比较上式两端可知$m_{34}\mathbf m_3=\mathbf r_3$，由于$\mathbf r_3$是正交单位矩阵的第三行，$\|\mathbf r_3\|=1$，那么有$m_{34}=\frac{1}{\|\mathbf m_3\|}$，因此利用$\mathbf M$矩阵可得摄像机内参
\begin{equation} \begin{aligned} u_0 & = \left(a_x \mathbf r_1^\top + u_0\mathbf r_3^\top\right)\mathbf r_3 = m_{34}^2\mathbf m_1^\top\mathbf m_3\\ v_0 & = \left(a_y \mathbf r_2^\top + v_0\mathbf r_3^\top\right)\mathbf r_3 = m_{34}^2\mathbf m_2^\top\mathbf m_3 \\ a_x & = m_{34}^2 \left\|\mathbf m_1\times\mathbf m_3\right\|\\ a_y & = m_{34}^2 \left\|\mathbf m_2\times\mathbf m_3\right\| \end{aligned}， \end{equation}
以及摄像机的外参
\begin{equation} \begin{aligned} \mathbf r_1 &= \frac{m_{34}}{a_x}\left(\mathbf m_1 - u_0\mathbf m_3\right)\\ \mathbf r_2 &= \frac{m_{34}}{a_y}\left(\mathbf m_2 - v_0\mathbf m_3\right)\\ \mathbf r_3 &= m_{34}\mathbf m_3\\ t_x &= \frac{m_{34}}{a_x}\left(m_{14} - u_0\right) \\ t_y &= \frac{m_{34}}{a_y}\left(m_{24} - v_0\right) \\ t_z &= m_{34} \end{aligned}。 \end{equation}
$\mathbf M$矩阵由4个摄像机内部参数及$\mathbf R$和$\mathbf t$确定。由于$\mathbf R$是正交单位阵，$\mathbf R$和$\mathbf t$的独立变量数为6，因此$\mathbf M$由10个独立变量的矩阵确定。但是$3\times 4$的矩阵$\mathbf M$由11个参数确定，可见这11参数并非相互独立，存在约束关系。在求解的时候并未参数间的约束关系，数据有误差时计算结果有误差，且误差在各参数间的分配也没有按约束关系考虑。实验表明，用上述方法分解出的内外参数有较大的误差。Faugeras等[4]在用\eqref{eq:ls-solution}求解$\mathbf M$时加上了约束条件$\|\mathbf m_3\|=1$，不需要解非线性方程，而且提高了计算精度。
改进的Faugeras标定法[3, P. 44]
像素畸变问题[5, P. 62]
由于制造工艺的限制，通常摄像机离散化后的像素不是矩形而是平行四边形，一边平行于$u$轴，另一边与$u$轴夹角为$\theta$，边长为分别为$\mathrm{d}x$和$\mathrm{d}y$，那么公式\eqref{eq:pixel-vs-image-coordination}可改写为
\begin{equation} \begin{bmatrix}u\\ v\\ 1\end{bmatrix} = \begin{bmatrix} {1\over \mathrm{d}x} & -{\cot\theta\over \mathrm{d}y} & u_0 \\ 0 & {\sin\theta\over \mathrm{d}y} & v_0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}x\\ y\\ 1\end{bmatrix}。 \label{eq:pixel-vs-image-coordination2} \end{equation}
Tsai标定法
实验表明，线性模型不能准确描述成像的几何关系，尤其在使用广角镜头时，远离图像中心处会有较大的畸变。描述非线性畸变可用如下公式[2, P. 59]
\begin{equation} \begin{aligned} \bar{x} &= x + \delta_x(x,y)\\ \bar{y} &= y + \delta_y(x,y) \end{aligned}， \label{eq:image-mapping-distortion} \end{equation}
其中，$(\bar x, \bar y)$是小孔线性模型计算出的图像坐标理想值，$(x, y)$是实际图像坐标，$\delta_x$和$\delta_y$是非线性畸变，它与图像点在图中位置相关，可用公式表示为
\begin{equation} \begin{aligned} \delta_x(x, y) &= k_1x\left(x^2+y^2\right)+\left(p_1\left(3x^2+y^2\right)+2p_2xy\right)+s_1\left(x^2+y^2\right)\\ \delta_y(x, y) &= k_2y\left(x^2+y^2\right)+\left(p_2\left(3x^2+y^2\right)+2p_1xy\right)+s_2\left(x^2+y^2\right) \end{aligned} \end{equation}
上式第一项称为径向畸变，第二项称为离心畸变（decentering），第三项称为薄棱镜畸变（thin prism），$k_1,k_2,p_1,p_2,s_1,s_2$称为非线性畸变参数。通常第一项径向畸变就能足够描述非线性畸变。Tsai指出[6]，考虑非线性畸变需要使用非线性优化算法，引入过多的非线性参数往往不能提高精度，反而引起解的不稳定性。但也有研究表明[7][8]，引入二三项在使用广角镜头时能提高精度。若只考虑径向畸变，\eqref{eq:image-mapping-distortion}可写为
\begin{equation} \begin{aligned} \bar{x} &= x\left(1+k_1r^2\right)\\ \bar{y} &= y\left(1+k_2r^2\right) \end{aligned}， \end{equation}
其中$r^2=x^2+y^2$，图像边缘处的畸变较大。
Tsai提出的方法[6][9]很好的解决了只存在径向畸变时的标定问题。Tsai标定法是基于径向排列约束（RAC，radial alignment constraint）的两步标定法。该方法的第一步是利用最小二乘法得到外部参数；第二步求解内部参数，如果摄像机无径向畸变，可由线性方程组求解，若存在径向畸变，可结合非线性优化得到全部参数。
张正友平面标定法
参考资料
- [1]Y. I. Abdel-Aziz and H. M. Karara, “Direct Linear Transformation from Comparator Coordinates into Object Space Coordinates in Close-Range Photogrammetry,” Photogrammetric Engineering & Remote Sensing, vol. 81, no. 2, pp. 103–107, 2015.
- [2]马颂德 and 张正友, 计算机视觉：计算理论与算法基础. 科学出版社, 1998.
- [3]徐德, 谭民, and 李原, 机器人视觉测量与控制, 2nd ed. 北京: 国防工业出版社, 2011.
- [4]O. D. Faugeras and G. Toscani, “The calibration problem for stereo,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 1986, vol. 86, pp. 15–20.
- [5]吴福朝, 计算机视觉中的数学方法. 科学出版社, 2008.
- [6]R. Y. Tsai, “An efficient and accurate camera calibration technique for 3D machine vision,” in Proc. IEEE Conf. on Computer Vision and Pattern Recognition, 1986.
- [7]O. D. Faugeras and G. Toscani, “Camera calibration for 3D computer vision,” in Proceedings of International Workshop on Machine Vision and Machine Intelligence, Tokyo, Japan, 1987.
- [8]J. Weng, P. Cohen, and M. Herniou, “Calibration of stereo cameras using a non-linear distortion model [CCD sensory],” in Pattern Recognition, 1990. Proceedings., 10th International Conference on, 1990, vol. 1, pp. 246–253.
- [9]R. Y. Tsai, “A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf TV cameras and lenses,” Robotics and Automation, IEEE Journal of, vol. 3, no. 4, pp. 323–344, 1987.
