相机内参标定

在网上找了一圈内参标定的博客，发现大多数都讲得太烂了，有的浅尝辄止，只讲了内外参是什么东西，有的更是离谱，写的人都没搞懂，不知道他们为什么要制造噪声。还是老老实实看论文吧，下面是1999年的正友标定法，论文名为：A Flexible New Technique for Camera Calibration。计算机行业多年轻，你看这么基础的东西才距离我们20年，所以不要灰心，还是能够学得会学得完的。

符号定义

先来一些符号定义，不要害怕，都是很简单的。

2D点写为$m=[u,v]^T$，3D点写为$M=[X,Y,Z]^T$，齐次坐标就是在上面加波浪号,也就是${\widetilde{\mathbf{m}}}=[u,v,1]^{T}$,${\widetilde{\mathbf{M}}}=[X,Y,Z,1]^{T}$。投影过程可写作：

$$ s\widetilde{\bf{m}}=\mathbf{A}\,[\mathbf{R}\ \ \ \ t]\widetilde{\bf{M}} $$

$s$是尺度因子，$R,t$是外参，$A$就是内参，这里投影公式和现在习惯的写法稍微有点不一样，但是原理是一致的，我还是尊重原文，方便后面的推导。$A$就是内参了，内参就是一个3x3的矩阵，如下所示：

$$ A=\begin{bmatrix}\alpha & \gamma & u_0 \\0 & \beta & v_0 \\0 & 0 & 1\end{bmatrix} $$

$(u_0,v_0)$是光心在像素坐标系上的投影，$\alpha,\beta$分别是$u,v$轴上的缩放因子，$\gamma$描述两个轴之间的倾斜，一般为0。后面如果出现了$A^{-T}$,它代表了$(A^{-1})^{T}$或者$(A^T)^{-1}$。

化简一下

待会儿我们用棋盘格做标定，所以点都在一个平面上，我们就假设一下，这个平面上所有点Z值为0，这是没问题的，因为世界坐标系本来就是任意设的。

$$ {S}\left[{\begin{array}{l}{u}\\ {v}\\ {1}\end{array}}\right]=\mathbf{A}\left[\mathbf{r}_{1}\quad\mathbf{r}_{2}\quad\mathbf{r}_{3}\quad\mathbf{t}\right]\left[{\begin{array}{l}{X}\\ {Y}\\ {0}\\ {1}\end{array}}\right] $$

$$ {S}\left[{\begin{array}{l}{u}\\ {v}\\ {1}\end{array}}\right]=\mathbf{A}\left[\mathbf{r}_{1}\quad\mathbf{r}_{2}\quad\mathbf{t}\right]\left[{\begin{array}{l}{X}\\ {Y}\\ {1}\end{array}}\right] $$

还套用之前的定义，只是这时候$M=[X,Y]^T$，${\widetilde{\mathbf{M}}}=[X,Y,1]^{T}$，有

$$ s\widetilde{\bf{m}}=H\widetilde{\bf{M}},\ \ \ \ H=\mathbf{A}\,\left[\mathbf{r}_{1}\quad\mathbf{r}_{2}\quad\mathbf{t}\right] $$

这么一来$H$就是一个3x3的矩阵,写开是下面这样，至于为什么多了一个$\lambda$，这是齐次矩阵的性质，整体乘以一个常数不影响，因为(1，2，1)和（2，4，2）代表了同一个点。有个前提知识是$H$是物理平面和成像平面间的单应性变换矩阵，这个是可以算出来的，利用平面上的点和成像平面点的匹配关系就能解，我博客有写单应性变换，不明白可以去瞅瞅。

$$ \left[\mathbf{h}_{1}\quad\mathbf{h}_{2}\quad\mathbf{h_3}\right]=\lambda\mathbf{A}\,\left[\mathbf{r}_{1}\quad\mathbf{r}_{2}\quad\mathbf{t}\right] $$

我们可以稍微化简一下上面的算式，可以看到如果内参$A$得到了，外参也能得到，$\lambda$反正是一个任意的数，我们这里重新定义$\lambda$，它其实是原来的$1/ \lambda$：

$$ r_1= \lambda A^{-1} h_1 $$

$$ r_2= \lambda A^{-1} h_2 $$

$$ r_3=r1 \times r2 $$

$$ t=\lambda A^{-1}h_3 $$

因为x轴，y轴正交，所以还有:

$$ {r_1}^Tr_2=0, \ \ \ {r_2}^Tr_1=0 $$

用代入法联合起来就是：

$$ \mathbf{h}_{1}^{T}\mathbf{A}^{-T}\mathbf{A}^{-1}\mathbf{h}_{2}=0 $$

$$ {\bf h}_{1}^{T}{\bf A}^{-T}{\bf A}^{-1}{\bf h}_{1}={\bf h}_{2}^{T}{\bf A}^{-T}{\bf A}^{-1}{\bf h}_{2} $$

这就是我们用来解内参的公式，当给定一个单应性矩阵时，我们可以得到两个基本的约束条件来限制内参。这是因为单应性矩阵有8个自由度，而外参有6个自由度，其中3个用于描述旋转，3个用于描述平移，这里多提一句，为什么旋转矩阵有9个参数，但自由度是3，那是旋转矩阵可以用三个欧拉角表示，所以9个参数间不是独立的。总而言之，一个单应性矩阵可以提供两个约束用来解内参。

解析法

内参是有解析解的，总结一下，步骤如下：

1.收集单应性矩阵：首先，你需要至少两个（通常是三个或更多）单应性变换矩阵。每个单应性变换矩阵对应于世界坐标系中的一个平面（棋盘格）在不同视角下的图像。为什么至少是两个呢？因为前面提到的，一个公式能给出两个约束也就可以解两个未知数，内参矩阵中的$\gamma$一般认为是0，所以内参矩阵有4个未知数；

2.构建方程：对于每个单应性变换矩阵$H_i$，可以构建以下方程，其中，$h_{i1}$和$h_{i2}$分别是单应性矩阵$H_i$的第一列和第二列，$A$是内参矩阵，$B = A^{-T} A^{-1}$：

$$ h_{i1}^T B h_{i2} = 0 ,\ \ \ h_{i1}^T B h_{i1} = h_{i2}^T B h_{i2} $$

3.线性求解：将上述方程组合起来，形成一个线性系统，可以通过线性求解方法（如奇异值分解）来求解B。

4.提取内参：一旦得到B，就可以通过对B进行分解来得到内参矩阵A。这通常涉及到对B进行乔里斯基（Cholesky）分解

但是呢，由于检测误差等各方面原因，算出来的内参矩阵$A$可能是不符合旋转矩阵的要求。

极大似然法

那我们就把误差考虑进去，假设我们得到了n张模型平面的图像，每个模型平面上有m个点。同时假设图像点受到独立同分布的噪声干扰，那么最大似然估计可以通过最小化以下函数获得：

$$ \begin{array}{l}{{\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}\big||{\bf{m}}_{ij}-\mathrm{\hat{m}}(A,\,{\bf{R}}_{i},\,t_{i},\,{\bf{M}}_{j})||^{2}}}\end{array} $$

这里

$$ \mathrm{\hat{m}}(A,\,{\bf{R}}_{i},\,t_{i},\,{\bf{M}}_{j}) $$

是$M_j$点在第$i$张图片上的投影。

重点来了我们可以通过前面的解析法得到初始估计的内外参，然后用这个初始内外参对所有点执行2D重投影，看看总体误差是多少，这里用的是非线性优化技术（如Levenberg-Marquardt算法）来精细调整内外参，以便最小化重投影误差。

畸变

上面其实已经讲完了，但是还不够完美，因为相机毕竟不是理想的针孔模型，不是投影矩阵能够描述的，比如通常的还有径向畸变，就是通常看到的鱼眼摄像头和运动相机那样，越远离中心越扭曲。假设$u,v$是无畸变的像素坐标，$\ddot{u},\ddot{v}$是实际观测的有畸变的，$x,y$是理想的相机坐标系坐标，$\ddot{x},\ddot{y}$则是实际观测到的有畸变的相机坐标系下坐标。径向畸变公式如下：

$$ \ddot{x}=x+x[k_{1}(x^{2}+y^{2})+k_{2}(x^{2}+y^{2})^{2}] $$

$$ {\ddot{y}}=y+y[k_{1}(x^{2}+y^{2})+k_{2}(x^{2}+y^{2})^{2}] $$

这里的$k_1,k_2$就是径向畸变的系数，径向畸变的中心就是主点。通过假设$\gamma$为0，利用相机坐标系到像素坐标系的转换有如下关系：

$$ \ddot{u}=u_0+\alpha \ddot{x} $$

$$ \ddot{v}=v_0+\beta \ddot{y} $$

代入到畸变公式：

$$ \ddot{u}=u+(u-u_0)[k_{1}(x^{2}+y^{2})+k_{2}(x^{2}+y^{2})^{2}] $$

$$ {\ddot{v}}=v+(v-v_0)[k_{1}(x^{2}+y^{2})+k_{2}(x^{2}+y^{2})^{2}] $$

下面就有两张方法来算畸变：

通过交替估计径向畸变

由于预期径向畸变会很小，人们会期望通过简单地忽略畸变，先利用我们前面的方法来算内参，然后再估计$k_1$与$k_2$，这推导出来的是理想的$u,v$，利用上面得到的畸变公式，每个图像的每个点都能解两个方程，如下所示：

$$ \left[ \begin{array}{cc} (u-u_0)(x^2+y^2) & (u-u_0)(x^2+y^2)^2\\ (v-v_0)(x^2+y^2) & (v-v_0)(x^2+y^2)^2 \\ \end{array} \right]\left[ \begin{array}{cc} k_1 \\ k_2\end{array}\right]=\left[ \begin{array}{cc}\ddot{u}-u\\ \ddot{v}-v \end{array}\right] $$

给定n张图像中的m个点，我们可以将所有方程堆叠起来，总共得到2mn个方程,用矩阵表示也就是$Dk=d$,这里$k=[k_1,k_2]^T$,线性最小二乘解如下：

$$ k=(D^TD)^{-1}D^Td $$

一旦估计出$k_1$和$k_2$，就又可以通过前面的极大似然估计来优化内外参，优化了内外参，再来优化畸变参数，不断重复这样交替优化的过程。

完整的最大似然估计

实验发现上述交替技术的收敛速度很慢。那么就把畸变参数融合进重投影，利用极大似然估计进行一把梭：

$$ \begin{array}{l}{{\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}\big||{\bf{m}}_{ij}-\mathrm{\hat{m}}(A,k_1,k_2,\,{\bf{R}}_{i},\,t_{i},\,{\bf{M}}_{j})||^{2}}}\end{array} $$

实际流程呢就是先把正常的投影过程做完继续做畸变，初始的畸变参数可以简单设为0，也可以用线性最小二乘解的结果作为初始值，每次同时优化内外参和畸变。

总结

这里真的是讲完了，概括一下步骤吧：

打印一个图案并将其贴在一个平面上；
通过移动平面或相机，从不同的方向拍摄模型平面的几张图片；
在图像中检测特征点；
利于解析法来估计五个内参和所有外参；
通过解线性最小二乘问题来估计径向畸变的系数；
把上述内外参和畸变参数作为初始值，通过最小化重投影误差来迭代改进所有参数。

棒棒，终于写完了，如果你觉得我写的有问题没明白，欢迎来博客主页的各种社交账号与我联系。