相机标定——张氏标定法

前言
- 动机
- 为什么要进行相机标定
- 什么是张氏标定法
张氏标定法的原理
- 透镜成像原理
- 世界坐标系到相机坐标系的转换
- 相机坐标系到图像坐标系的转换
- 图像坐标系到像素坐标系
- 单应性矩阵
- 内参求解
- 外参求解
- Matlab实操

前言

动机

前段时间在整机械臂的手眼标定，也就是标定3D相机到机械臂末端坐标系之间的转化关系，需要用到Matlab的相机标定工具来得到相机与标定板之间的齐次变换矩阵（也就是相机标定中的外参）。Matlab标定工具箱使用起来确实方便，不需要知道其中的标定原理，但是我认为要使用好一个工具，对原理有一定的了解是必要的，故而花了点时间研究了一下Matlab标定工具箱使用的标定算法——张氏标定法。

为什么要进行相机标定

相机需要标定，总结起来主要有以下两点原因：

由于透镜制造精度以及组装工艺的偏差会引入畸变，导致原始图像失真，需要进行矫正

在进行一些视觉任务如三维重建、视觉抓取的时候，我们需要建立从世界坐标系到图像坐标系之间的映射关系。例如可以通过某一点P在图像坐标系中的位置推导出其在世界坐标系中的位置，这是进行后续工作的基础。

什么是张氏标定法

张氏标定法是张正友博士在1999年发表在国际顶级会议ICCV上的论文《Flexible Camera Calibration By Viewing a Plane From Unknown Orientations》中，提出的一种利用平面棋盘格进行相机标定的实用方法。其后2000的这篇《A flexible new technique for camera calibration》引用数更是达到了恐怖的13885次，江湖地位可见一斑。该方法介于摄影标定法和自标定法之间，既克服了摄影标定法需要的高精度三维标定物(贵，操作麻烦)的缺点，又解决了自标定法鲁棒性差的难题。标定过程仅需使用一个打印出来的棋盘格，并从不同方向拍摄几组图片即可（理论上3组就可以求解所有的未知数，不过一般采15~20组，用优化的方法保证求解精度较高），任何人都可以自己制作标定图案，不仅实用灵活方便，而且精度很高，鲁棒性好。因此很快被全世界广泛采用，极大的促进了三维计算机视觉从实验室走向真实世界的进程。（摘自张正友标定法-完整学习笔记-从原理到实战 - 知乎 (zhihu.com)）
张氏标定法相比于其他标定方法，在保证一定精度以及鲁棒性的情况下，操作简单、成本较低。但是有一点要注意，张氏标定法只考虑了径向畸变，没有考虑切向畸变。
Tips(本文为了方便（简单），不考虑相机的畸变)

张氏标定法的原理

张氏标定法的整体思路为先求出世界坐标系到像素坐标系的单应性矩阵（单应性矩阵其实就是一个图像中的像素点到另外一个图像中像素点的变换矩阵，双目相机系统中也有单应性矩阵，其是左右相机图像之间的变换矩阵），然后根据单应性矩阵得到内参矩阵，最后得到外参矩阵，不过这其中有个东西要考虑，就是尺度因子，我花了一段比较长的时间去思考这个。

透镜成像原理

我将主要用图来呈现，比码字快而且直观，透镜成像原理如下所示：

上图总的来说就是建立了从像素坐标系到世界坐标系之间的映射关系。其步骤为:先通过外参建立相机坐标系与世界坐标系之间的联系，然后通过内参与相机成像模型建立相机坐标系与图像坐标系，图像坐标系与像素坐标系之间的联系。接下来，我们对上图的步骤进行拆解与分析。

世界坐标系到相机坐标系的转换

上式子中，[x_c,y_c,z_c] ， [x_w,y_w,z_w] 分别为点在相机坐标系与世界坐标系下的位置（注意，这还没有投影变换，只是换了个参考坐标系，R为从世界坐标系到相机坐标系的旋转矩阵，T为平移向量，两者构成了4*4的齐次变换矩阵，这也就是标定过程中得到的外参，这个外参正常情况下可能用不上，但是在机械臂的手眼标定中这个外参是必须的参数之一。

相机坐标系到图像坐标系的转换

相机坐标系到图像坐标系的关系如下图所示：
从相机坐标系到图像坐标系经过了投影变换，图像坐标为二维坐标，没有Z轴分量，成像平面距光心的距离为焦距f,这边有个要注意的是，从相机坐标系到图像坐标系的变化会引入一个尺度因子，这个尺度因子和 Z_c 相关。

图像坐标系到像素坐标系

相机实际输出的图片是按像素进行排列的，我们最后得到的数据也是像素的位置，像素坐标系相当于对图像坐标系做了原点偏移与坐标值离散化。

其中，1/dx, 1/dy 分别为像素在 x , y 轴方向上的物理尺寸，u₀ , v₀ 为主点（图像原点）在像素平面的偏移。（其实正常来说的话，在上图3*3矩阵中的第一行第二列位置还会有一个 γ ，这个参数是 x, y 轴之间的主轴偏差，因为实际透镜，两个光轴可能不是垂直的）

单应性矩阵

接下来，我们综合上述从世界到相机、从相机到图像、从图像到像素的转化关系，可以得到如下式子：

我们对式子进行化简，我们可以将世界坐标系的 XOY 平面移动到标定板上，那么标定板上的点在世界坐标系上的 z 坐标就是0，显然，这只是换了一种世界坐标系的选取方式，只会影响外参，对内参没有任何影响。在此基础上，进行合并，我们可得：

其中，s 为尺度因子，其之后分别为内参矩阵与外参矩阵。这儿，我们就可以得到单应性矩阵 H 了,其表达式如下所示

可以看出来，单应性矩阵同时包含外参与内参，同时还含有尺度因子，求解的一般步骤是先求单应性矩阵，然后再求内外参，求解步骤如下：
我们先把 H 矩阵展开，如下：

我们假设像素坐标系与世界坐标系对应点对齐次坐标为(x’,y’,1)和(x,y,1)（z=0），那么有

展开之后，可以得到如下式子：

注意，上述两个式子为齐次方程，h_ij 同时除以一个倍数对结果没有影响，（这也可以看出来，单应性矩阵的自由度为8，这边给出单应性矩阵自由度为8的解释：如果给定一个单应H={h_ij}，给它的元素乘上同一个数a，得到的的单应aH和H作用相同，因为新单应无非把齐次点x1变成了齐次点ax1，都是一回事。因此我们可以把a换成1/h33，那么H就变成了只有8个自由元素的矩阵），注意，按上述方法求解相当于对单应矩阵做了归一化类似的操作，故而求解出来的矩阵与真正的H矩阵之间还差一个倍数关系，这个倍数关系与我们加约束相关，也就是说如果我们加约束h33=1,或者h32=1，倍数是不一样的。
由于单应矩阵有8个自由度（我们加一个约束，h₃₃=1,其他形式的约束也可以），我们需要同一个图像上的四个点便可求解出单应性矩阵，当然，点越多越精确，求解方法如下：

如果我们得到了N个点对，有：

可以采用最小二乘法等优化的方法求解，注意，这边求出来的 H 矩阵不可以直接用于世界坐标系与像素坐标系的对应，我们求出来的是一个做了归一化类似操作的矩阵。

内参求解

单应性矩阵如下所示：

接下来我们进行内参的求解，步骤如下：
按元素对应关系，可以得到如下式子：

这边我们得注意，r1与r2是什么，这两个是旋转向量，旋转向量有两个性质

两个旋转向量之间相互正交
旋转向量的模为1

我们可以根据这两个性质构建约束
约束1：

约束2：

为了方便，我们构造B，B的表达形式如下：

可以看出来，B其实是对称矩阵，我们只需要求出上三角的六个元素便可，我们把H矩阵的列向量h_i展开，有

那么，代入约束1，我们可得：

其实这就是把左边的式子展开，然后合并同类项就行，左边式子展开就是一个多项式，右边的展开也是一个多项式。
为了简化表达，我们令：

那么最后的约束条件可与转化成：

v₁₂,v₁₁,v₂₂可以通过前面已经计算好的单应矩阵得到，是已经知道的，那么我们只需要采集3张以上的图片，就可以求解B矩阵，进通过如下式子求解内参：

这一块我有一个疑惑，求内参需要考虑这个 λ 么，我感觉不需要考虑这个，它一直处于被消去的状态，上述公式中的λ 是不是就是1就行（望大神赐教）

外参求解

接下来，我们对外参进行求解，外参是相对于每一张标定图片而言的，不同的图片也就对应了不同的外参。

其中又由旋转矩阵性质有：

则可得：

那么，r₁,r₂r₃,t便都可以求出来了。

看论文以及其他blog的时候，主要就是内参求解的 λ 这一块我没太看懂，我不知道我上述的理解对不对，我感觉内参求解确实是不需要考虑尺度因子的

Matlab实操

这块不想写了，网上资料很多，请看链接：MATLAB–相机标定教程
有几个东西我解释一下吧,输入黑白格边长是为了计算像素点世界坐标用的，所以这个的准确性势必会影响到标定结果，原点的话，MATLAB会自动确定的，不过你得保证你标定的时候，所有原点都在一个点上（保证标定板长宽格子数不对称一般就不会出现这个问题），采图15~20 即可,标定板视野占图片的1/4-1/3。

期待能有机会与各位大佬交流学习，如有错误，欢迎指正！！！