Camera Calibration-526互联

Camera Calibration Algorithm

前言

首先要了解知道什么是相机的针孔模型，可参考文章小孔相机参数学习笔记。

先介绍相机标定。所谓相机标定，就是通过实验等方法将相机模型中的未知数给解算出来。例如小孔模型我们一般需要节算出内参矩阵 fx，fy，cx，cy 以及畸变参数。

为什么要标定出这些参数呢？
一个是因为每个镜头的畸变程度各不相同，通过相机标定可以校正这种镜头畸变矫正畸变，生成矫正后的图像；另一个是根据获得的图像重构三维场景，例如__ PNP __算法（openCV中solvePNP函数），就是通过相机的内参以及畸变解算出目标在相机下的位置，以及双目 3D 相机也需要相机模型参数进行计算等。

有很多标定方法，比如传统相机标定法、主动视觉相机标定方法、相机自标定法等。现在用的最多的是相机的自标定法，它的基础即张正友标定法。

“张正友标定法” 是张正友博士在 1999 年发表在国际顶级会议 ICCV 上的论文《Flexible Camera Calibration By Viewing a Plane From Unknown Orientations》中，提出的一种利用平面棋盘格进行相机标定的实用方法。该方法既克服了摄影标定法需要的高精度三维标定物的缺点，又解决了之前自标定法鲁棒性差的难题。

标定过程仅需使用一个打印出来的棋盘格，并从不同方向拍摄几组图片即可，任何人都可以自己制作标定图案，不仅实用灵活方便，而且精度很高，鲁棒性好。因此很快被全世界广泛采用。

张正友标定法计算流程

首先我们看一下张正友标定法使用 OpenCV 的计算流程：

准备标定图片，原理上三张就够，一般在多个角度采集 20 张左右。
提取标定板的关键点，并计算出标定板上关键点的实际相对位置，一般将标定板当做 XY 平面，Z 为 0，标定板第一个点为坐标原点。
相机标定，通过张正友标定法计算出内参外参以及畸变。
对标定结果进行评价，一般通过重投影的误差进行评价。
查看标定效果，利用标定结果对棋盘图进行矫正

从上边的过程可以看出，其实只有第三步是真正的解算过程，大致的方法如下：

首先用于标定的棋盘格是三维场景中的一个平面Π，其在成像平面的像是另一个平面?，知道了两个平面的对应点的坐标，就可以求解得到两个平面的单应矩阵?。其中，标定的棋盘格是特制的，其角点的坐标是已知的；图像中的角点，可以通过角点提取算法得到，这样就可以得到棋盘平面Π和图像平面?的单应矩阵?，即：

\[p=K[R|t]P \]

其中p是像素点坐标，P是标定板上的棋盘格坐标，K 是相机内参，为了不失一般性，可以在相机的内参矩阵上添加一个扭曲参数γ。即:

\[K=\left[\begin{array}{ccc}f_x&&\gamma&&c_x\\ 0&&f_y&&c_y\\ 0&&0&&1\end{array}\right] \]

这里假设H为单应矩阵，通常用于表示两个平面之间的投影变换，即p=HP。这样就可以得到下面的等式：

\[H=K[R|t] \]

我们通过对应的点对解得H后（即p=HP），则可以通过上面的等式得到相机的内参数K，以及外参旋转矩阵R和平移向量t。

下面具体推导具体的解法：

设棋盘格所在的平面为世界坐标系中Z=0 的平面，这样棋盘格的任一角点P的世界坐标为 (X,Y,0)，根据小孔相机模型（小孔成像中s表示镜头到被成像物体的距离，镜头当作一个点处理；f为相机的焦距，即镜头到感光元件的距离。那么根据相似三角形，即有p/f=P/s，也即sp=fP。这里f焦距即内参，是矩阵K的一部分，就有了下面的式子。）：

\[s\begin{pmatrix}u\\ v\\ 1\end{pmatrix}=K\begin{bmatrix}R&t\end{bmatrix}\begin{pmatrix}X\\ Y\\ 0\\ 1\end{pmatrix}=K[r_1\quad r_2\quad r_3\quad t]\begin{pmatrix}X\\ Y\\ 0\\ 1\end{pmatrix}=K[r_1\quad r_2\quad t]\begin{pmatrix}X\\ Y\\ 1\end{pmatrix} \]

另外这里需要在列向量(X Y 0)（Z为0）上增加一维“1”，以添加[R t]坐标系转换矩阵中t平移的部分。这里肯定单靠旋转不足以变换坐标系的，需要平移向量t的信息。那么等号左侧列向量(u v)也需要加“1”变成(u v 1)，这样左边是一个3X1的列向量，右边是K-3X3乘[R t]-3X4再乘(X Y 0 1)-4X1，即可得到3X1列向量。

再根据单应性原则：

\[H=\begin{bmatrix}h_1&h_2&h_3\end{bmatrix}=\begin{array}{c c c}{{=\lambda K[r_{1}}}&{{r_{2}}}&{{t]}}\end{array} \]

根据上面两式则有：

\[\begin{gathered} r_{1}=\lambda K^{-1}h_{1} \\ r_{2}=\lambda K^{-1}h_{2} \\ t=\lambda K^{-1}h_{3} \end{gathered} \]

将旋转矩阵R的各个列向量和平移向量t使用H的列向量表示：

又因为，R是旋转矩阵，则其是正交矩阵，也就是其任意两个列向量的内积为 0，列向量的模为 1，则有：

即：

\[\left\{\begin{array}{l}{{h_{1}^{T}K^{-T}K^{-1}h_{2}=0}}\\ {{h_{1}^{T}K^{-T}K^{-1}h_{1}=h_{2}^{T}K^{-T}K^{-1}h_{2}=1}}\end{array}\right. \]

$\begin{aligned}&r_{1}^{T}r_{2}=0\&\\&|r_{1}|=|r_2|=1\end{aligned}$那么对于一幅棋盘标定版的图像（一个单应矩阵）可以获得两个对内参数的约束等式。(H由对应点可得，求K即可。）

我们令：

\[B=K^{-T}K^{-1}=\begin{bmatrix}B_{11}&B_{12}&B_{13}\\ B_{21}&B_{22}&B_{23}\\ B_{31}&B_{32}&B_{33}\end{bmatrix} \]

矩阵B是一个对称矩阵，其未知量只有 6 个，将 6 个未知量写为向量的形式：

\[b=\left[\begin{array}{c}B_{11},B_{12},B_{22},B_{13},B_{23},B_{33}\end{array}\right] \]

则有：

\[h_i K^{-T}K^{-1}h_j=h_i B h_j=v_{ij}^T b \]

其中：

\[v_{ij}=\begin{bmatrix}h_{i1}h_{j1}&h_{i1}h_{j2}+h_{i2}h_{j1}&h_{i2}h_{j2}&h_{i3}h_{j1}+h_{i1}h_{j3}&h_{i3}h_{j2}+h_{i2}h_{j3}&h_{i3}h_{j3}\end{bmatrix}^T \]

则约束等式有：

\[\left\{\begin{array}{l}{{h_{1}^{T}K^{-T}K^{-1}h_{2}=0}}\\ {{h_{1}^{T}K^{-T}K^{-1}h_{1}=h_{2}^{T}K^{-T}K^{-1}h_{2}=1}}\end{array}\right.\Rightarrow\left\{\begin{array}{l}{{v_{22}^{T}b=0}}\\ {{v_{11}b=v_{12}b}}\end{array}\right. \]

写成矩阵的形式有：

\[\left[\begin{array}{c}v_{12}^T\\ v_{11}-v_{22}\end{array}\right]b=0 \]

假如有n幅图像，则:

\[\mathrm{~V~}b =0 \]

其中，V是一个 2n×6 的矩阵，b是一个 6 维向量，所以

当n≥3，可以得到b的唯一解;
当n=2，则可以假设扭曲参数γ=0 作为额外的约束条件
当n=1，则值能计算两个相机的内参数

对于方程Vb=0 可以使用 SVD 求得其最小二乘解。对V?V进行 SVD 分解，其最小特征值对应的特征向量就是Vb=0 的最小二乘解，从而求得矩阵B。由于这里得到的B的估计值是在相差一个常量因子下得到的，所以有：

\[B=\lambda A^{-T}A \]

所以则有：

$\begin{cases}c_x=\gamma c_y/\beta-B_{13}\alpha^2/\lambda\\c_y=(B_{12}B_{11}-B_{11}B_{23})/(B_{11}B_{22}-B_{12}^2)\\\alpha=\sqrt{\lambda B_{11}/B_{11}}\\\beta=\sqrt{\lambda B_{11}/B_{11}B_{22}-B_{12}^2)}\\\gamma=-B_{12}\alpha^2/\lambda\\\lambda=B_{33}-[B_{13}^2+c_3(B_{12}B_{13}-B_{11}B_{23})]/B_{11}\end{cases}$

其中 $fx = α_(1/γ),fy = β_(1/γ) $。

为了进一步增加标定结果的可靠性，可以使用最大似然估计来优化上面估计得到的结果。

假设同一相机从n个不同的角度的得到了n幅标定板的图像，每幅图像上有m个像点。Mi?表示第i幅图像上第j个像点对应的标定板上的三维点，则：

\[\hat{m}(K,R_i,t_i,M_{ij})=K[R\quad t]M_{ij} \]

m̂ (K,Ri,ti,Mij) 表示Mij的像点。其中，Ri,ti表示第i幅图像对应相机的旋转矩阵和平移向量，K是相机的内参数。则像点mij的概率密度函数是：

\[f(m_{ij})=\frac{1}{\sqrt{2\pi}}e^{\frac{-(\hat{m}(K,R_{i},t_{i},M_{ij})-m_{ij})^{2}}{\sigma^{2}}} \]

构造似然函数：

\[L(A,R_{i},t_{i},M_{i j})=\prod_{i=1,j=1}^{n,m}f(m_{i j})={\frac{1}{\sqrt{2\pi}}}e^{{\frac{-\Sigma_{i=1}^{n}\sum_{j=1}^{m}(\hat m( K,R_{i}, t_{i},M_{i j})-m_{i j})^{2}}{\sigma^{2}}}} \]

为了能够让L取得最大值，需要最小化下面的值

$\sum\limits_{i=1}^{n}\|\hat{m}(K,R_i,t_i,M_{ij})-m_{ij}\|^2$

问题变成了一个非线性优化问题，利用上面得到的解作为初始值，迭代得到最优解。这个过程就是在减少重投影误差的过程。

至此，通过张正友标定法，我们获得了相机的内参以及外参，但是畸变没有获得。张正友标定法只关注了影响较大的径向畸变。畸变的解算有点类似内参解算，暂时先不列举。