尊重别人的劳动成果就是对自己的尊重——声明至上：转载来源https://zhuanlan.zhihu.com/p/35223115
从零开始学习「张氏相机标定法」（一）相机成像模型
先来简单介绍一下我们的主角：张正友博士。他是世界著名的计算机视觉和多媒体技术的专家，ACM Fellow，IEEE Fellow。现任微软研究院视觉技术组高级研究员。他在立体视觉、三维重建、运动分析、图像配准、摄像机标定等方面都有开创性的贡献。
「张氏标定法」是张正友博士在1999年发表在国际顶级会议ICCV上的论文《Flexible Camera Calibration By Viewing a Plane From Unknown Orientations》中，提出的一种利用平面棋盘格进行相机标定的实用方法。该方法介于摄影标定法和自标定法之间，既克服了摄影标定法需要的高精度三维标定物的缺点，又解决了自标定法鲁棒性差的难题。标定过程仅需使用一个打印出来的棋盘格，并从不同方向拍摄几组图片即可，任何人都可以自己制作标定图案，不仅实用灵活方便，而且精度很高，鲁棒性好。因此很快被全世界广泛采用，极大的促进了三维计算机视觉从实验室走向真实世界的进程。
在介绍「张氏标定法」之前，我们先来搞清楚一些基本问题。
问题一：为什么要进行相机标定？
相机标定的目的是：建立相机成像几何模型并矫正透镜畸变。这句话有点拗口，下面分别对其中两个关键部分进行解释。
建立相机成像几何模型：计算机视觉的首要任务就是要通过拍摄到的图像信息获取到物体在真实三维世界里相对应的信息，于是，建立物体从三维世界映射到相机成像平面这一过程中的几何模型就显得尤为重要，而这一过程最关键的部分就是要得到相机的内参和外参（后续文有具体解释）。
矫正透镜畸变：我们最开始接触到的成像方面的知识应该是有关小孔成像的，但是由于这种成像方式只有小孔部分能透过光线就会导致物体的成像亮度很低，于是聪明的人类发明了透镜。虽然亮度问题解决了，但是新的问题又来了：由于透镜的制造工艺，会使成像产生多种形式的畸变，于是为了去除畸变（使成像后的图像与真实世界的景象保持一致），人们计算并利用畸变系数来矫正这种像差。虽然理论上可以设计出不产生畸变的透镜，但其制造工艺相对于球面透镜会复杂很多，所以相对于复杂且高成本的制造工艺，人们更喜欢用数学来解决问题。
相机成像模型
前面已经说过，相机标定的目的之一是为了建立物体从三维世界到成像平面上各坐标点的对应关系，所以首先我们需要定义这样几个坐标系来为整个过程做好铺垫：
世界坐标系(world coordinate system)：用户定义的三维世界的坐标系，为了描述目标物在真实世界里的位置而被引入。单位为m。
相机坐标系(camera coordinate system)：在相机上建立的坐标系，为了从相机的角度描述物体位置而定义，作为沟通世界坐标系和图像/像素坐标系的中间一环。单位为m。
图像坐标系(image coordinate system)：为了描述成像过程中物体从相机坐标系到图像坐标系的投影透射关系而引入，方便进一步得到像素坐标系下的坐标。 单位为m。
像素坐标系(pixel coordinate system)：为了描述物体成像后的像点在数字图像上（相片）的坐标而引入，是我们真正从相机内读取到的信息所在的坐标系。单位为个（像素数目）。
一下子定义出来四个坐标系可能有点晕，下图可以更清晰地表达这四个坐标系之间的关系：
世界坐标系：Xw、Yw、Zw。相机坐标系： Xc、Yc、Zc。图像坐标系：x、y。像素坐标系：u、v。
其中，相机坐标系的 轴与光轴重合，且垂直于图像坐标系平面并通过图像坐标系的原点，相机坐标系与图像坐标系之间的距离为焦距f（也即图像坐标系原点与焦点重合）。像素坐标系平面u-v和图像坐标系平面x-y重合，但像素坐标系原点位于图中左上角（之所以这么定义，目的是从存储信息的首地址开始读写）。
在这里我们先引入「棋盘」的概念：
棋盘是一块由黑白方块间隔组成的标定板，我们用它来作为相机标定的标定物（从真实世界映射到数字图像内的对象）。之所以我们用棋盘作为标定物是因为平面棋盘模式更容易处理（相对于复杂的三维物体），但与此同时，二维物体相对于三维物体会缺少一部分信息，于是我们会多次改变棋盘的方位来捕捉图像，以求获得更丰富的坐标信息。如下图所示，是相机在不同方向下拍摄的同一个棋盘图像。
下面将依次对刚体进行一系列变换，使之从世界坐标系进行仿射变换、投影透射，最终得到像素坐标系下的离散图像点，过程中会逐步引入各参数矩阵。
1、从世界坐标系到相机坐标系
刚体从世界坐标系转换到相机坐标系的过程，可以通过旋转和平移来得到，我们将其变换矩阵由一个旋转矩阵和平移向量组合成的齐次坐标矩阵（为什么要引入齐次坐标可见后续文章）来表示：
其中，R为旋转矩阵，t为平移向量，因为假定在世界坐标系中物点所在平面过世界坐标系原点且与Zw轴垂直（也即棋盘平面与Xw-Yw平面重合，目的在于方便后续计算），所以zw=0，可直接转换成式1的形式。其中变换矩阵
即为前文提到的外参矩阵，之所称之为外参矩阵可以理解为只与相机外部参数有关，且外参矩阵随刚体位置的变化而变化。
下图表示了用R，t将上述世界坐标系转换到相机坐标系的过程。
2、从相机坐标系到理想图像坐标系（不考虑畸变）
这一过程进行了从三维坐标到二维坐标的转换，也即投影透视过程（用中心投影法将物体投射到投影面上，从而获得的一种较为接近视觉效果的单面投影图，也就是使我们人眼看到景物近大远小的一种成像方式）。我们还是拿针孔成像来说明（除了成像亮度低外，成像效果和透镜成像是一样的，但是光路更简单）。
成像过程如图二所示：针孔面（相机坐标系）在图像平面（图像坐标系）和物点平面（棋盘平面）之间，所成图像为倒立实像。
但是为了在数学上更方便描述，我们将相机坐标系和图像坐标系位置对调，变成图三所示的布置方式（没有实际的物理意义，只是方便计算）：
此时，假设相机坐标系中有一点M，则在理想图像坐标系下（无畸变）的成像点P的坐标为（可由相似三角形原则得出）：
将上式化为齐次坐标表示形式为：
3、从理想图像坐标系到实际图像坐标系（考虑畸变）
透镜的畸变主要分为径向畸变和切向畸变，还有薄透镜畸变等等，但都没有径向和切向畸变影响显著，所以我们在这里只考虑径向和切向畸变。
径向畸变是由于透镜形状的制造工艺导致。且越向透镜边缘移动径向畸变越严重。下图所示是径向畸变的两种类型：桶形畸变和枕形畸变。
实际情况中我们常用r=0处的泰勒级数展开的前几项来近似描述径向畸变。矫正径向畸变前后的坐标关系为：
由此可知对于径向畸变，我们有3个畸变参数需要求解。
切向畸变是由于透镜和CMOS或者CCD的安装位置误差导致。因此，如果存在切向畸变，一个矩形被投影到成像平面上时，很可能会变成一个梯形。切向畸变需要两个额外的畸变参数来描述，矫正前后的坐标关系为：
由此可知对于切向畸变，我们有2个畸变参数需要求解。
综上，我们一共需要5个畸变参数（k1、k2、k3、p1和p2 ）来描述透镜畸变。
4、从实际图像坐标系到像素坐标系
由于定义的像素坐标系原点与图像坐标系原点不重合，假设图像坐标系原点在像素坐标系下的坐标为（u0，v0），每个像素点在图像坐标系x轴、y轴方向的尺寸为：dx、dy，且像点在实际图像坐标系下的坐标为（xc，yc），于是可得到像点在像素坐标系下的坐标为：
化为齐次坐标表示形式可得：
公式2中(xp, yp)与公式5中(xc, yc)相同，都是图像坐标系下的坐标。
若暂不考虑透镜畸变，则将式2与式5的转换矩阵相乘即为内参矩阵M：
之所以称之为内参矩阵可以理解为矩阵内各值只与相机内部参数有关，且不随物体位置变化而变化。
最后用一幅图来总结从世界坐标系到像素坐标系（不考虑畸变）的转换关系：
从零开始学习「张氏相机标定法」（二）单应矩阵
前面文章《从零开始学习「张氏相机标定法」（一）成像几何模型》中我们已经得到了像素坐标系和世界坐标系下的坐标映射关系：
其中，u、v表示像素坐标系中的坐标，s表示尺度因子，fx、fy、u0、v0、γ（由于制造误差产生的两个坐标轴偏斜参数，通常很小）表示5个相机内参，R,t表示相机外参，Xw、Yw、Zw（假设标定棋盘位于世界坐标系中Zw=0的平面）表示世界坐标系中的坐标。
单应性概念的引出
我们在这里引入一个新的概念：单应性（Homography）变换。可以简单的理解为它用来描述物体在世界坐标系和像素坐标系之间的位置映射关系。对应的变换矩阵称为单应性矩阵。在上述式子中，单应性矩阵定义为：
其中，M是内参矩阵
从单应矩阵定义式子来看，它同时包含了相机内参和外参。在进一步介绍相机标定知识之前，我们重点来了解一下单应性，这有助于深入理解相机标定。因为在计算机视觉领域，单应性是一个非常重要的概念。
为了不让读者一上来就淹没在公式的汪洋大海中失去兴趣，我们颠倒一下顺序，先来看看单应性到底有什么用，然后再介绍单应矩阵的估计方法。
单应性在计算机视觉中的应用
单应性在计算机视觉领域是一个非常重要的概念，它在图像校正、图像拼接、相机位姿估计、视觉SLAM等领域有非常重要的作用。
1、图像校正
用单应矩阵进行图像矫正的例子如下图所示，最少需要四个对应点对（后面会给出原因）就可以实现。
2、视角变换
单应矩阵用于视角变换的例子如下图所示，可以方便地将左边普通视图转换为右图的鸟瞰图。
3、图像拼接
既然单应矩阵可以进行视角转换，那我们把不同角度拍摄的图像都转换到同样的视角下，就可以实现图像拼接了。如下图所示，通过单应矩阵H可以将image1和image2都变换到同一个平面。
单应矩阵用于图像拼接的例子如下所示。
4、增强现实（AR）
平面二维标记图案（marker）经常用来做AR展示。根据marker不同视角下的图像可以方便的得到虚拟物体的位置姿态并进行显示，如下图所示。
如何估计单应矩阵？
了解了上述单应性的部分应用后，我们就有很大的动力来学习单应矩阵的推导和计算了。首先，我们假设两张图像中的对应点对齐次坐标为(x',y',1)和(x,y,1)，单应矩阵H定义为：
则有：
矩阵展开后有3个等式，将第3个等式代入前两个等式中可得：
也就是说，一个点对对应两个等式。在此插入一个讨论：单应矩阵H有几个自由度？
或许有人会说，9个啊，H矩阵不是9个参数吗？从h11到h33总共9个。真的是这样吗？实际上并不是，因为这里使用的是齐次坐标系，也就是说可以进行任意尺度的缩放。比如我们把hij乘以任意一个非零常数k并不改变等式结果：
所以实际上单应矩阵H只有8个自由度。8自由度下H计算过程有两种方法。
第一种方法：直接设置 h33=1，那么上述等式变为：
第二种方法：将H添加约束条件，将H矩阵模变为1，如下：
以第2种方法（用第1种也类似）为例继续推导，我们将如下等式（包含||H||=1约束）:
乘以分母展开，得到：
整理，得到：
假如我们得到了两幅图片中对应的N个点对（特征点匹配对），那么可以得到如下线性方程组：
写成矩阵形式：
由于单应矩阵H包含了||H||=1约束，因此根据上图的线性方程组，8自由度的H我们至少需要4对对应的点才能计算出单应矩阵。这也回答了前面图像校正中提到的为何至少需要4个点对的根本原因。
但是，以上只是理论推导，在真实的应用场景中，我们计算的点对中都会包含噪声。比如点的位置偏差几个像素，甚至出现特征点对误匹配的现象，如果只使用4个点对来计算单应矩阵，那会出现很大的误差。因此，为了使得计算更精确，一般都会使用远大于4个点对来计算单应矩阵。另外上述方程组采用直接线性解法通常很难得到最优解，所以实际使用中一般会用其他优化方法，如奇异值分解、Levenberg-Marquarat（LM）算法（后续文章会介绍）等进行求解。
如何根据标定图得到单应矩阵？
经过前面一系列的介绍，我们应该大致明白如何根据打印的棋盘标定图和拍摄的照片来计算单应矩阵H。我们来总结一下大致过程。
1、打印一张棋盘格标定图纸，将其贴在平面物体的表面。
2、拍摄一组不同方向棋盘格的图片，可以通过移动相机来实现，也可以移动标定图片来实现。
3、对于每张拍摄的棋盘图片，检测图片中所有棋盘格的特征点（角点，也就是下图中黑白棋盘交叉点，中间品红色的圆圈内就是一个角点）。我们定义打印的棋盘图纸位于世界坐标系Zw=0的平面上，世界坐标系的原点位于棋盘图纸的固定一角（比如下图中黄色点）。像素坐标系原点位于图片左上角。
4、因为棋盘标定图纸中所有角点的空间坐标是已知的，这些角点对应在拍摄的标定图片中的角点的像素坐标也是已知的，如果我们得到这样的N>=4个匹配点对（越多计算结果越鲁棒），就可以根据LM等优化方法得到其单应矩阵H。当然计算单应矩阵一般不需要自己写函数实现，OpenCV中就有现成的函数可以调用，对应的c++函数是:
Mat findHomography(InputArray srcPoints, InputArray dstPoints, int method=0, double ransacReprojThreshold=3, OutputArray mask=noArray() )
从函数定义来看，只要输入匹配点对，指定具体计算方法即可输出结果。
至此，我们已经搞清楚单应矩阵的概念、推导和应用了，下一篇文章我们继续学习张氏相机标定法。
从零开始学习「张氏相机标定法」（三）推导求解
前面文章《从零开始学习「张氏相机标定法」（一）成像几何模型》中我们学习了相机成像几何模型，知道了如何将世界坐标系中的三维坐标和像素坐标系中的二维坐标联系起来，根据《从零开始学习「张氏相机标定法」（二）单应矩阵》，我们根据标定棋盘图纸及其对应的照片已经可以得到单应矩阵H了。如下所示：
下一步如何求相机内外参数呢？
我们知道H是内参矩阵和外参矩阵的混合体，而我们想要最终分别获得内参和外参。所以需要想个办法，先把内参求出来（先求内参是因为更容易），得到内参后，外参也就随之解出了。
我们先不考虑镜头畸变，来看看如何求解内参和外参。求解思路是利用旋转向量的约束关系，以下是具体推导，建议自己演算一遍，加深理解。
为了利用旋转向量之间的约束关系，我们先将单应性矩阵H化为3个列向量，即H=[h1 h2 h3]，则有
按元素对应关系可得：
因为旋转向量在构造中是相互正交的，即r1和r2相互正交，由此我们就可以利用“正交”的两个含义，得出每个单应矩阵提供的两个约束条件：
约束条件1：旋转向量点积为0（两垂直平面上的旋转向量互相垂直），即：
约束条件2：旋转向量长度相等（旋转不改变尺度），即：
所以一个单应性矩阵H可以提供上述两个约束条件。那么如何利用上述两个约束条件求解内参或者外参呢？我们一步一步来看，由前面可知内参矩阵M：
记：
我们看到B为对称矩阵，真正有用的元素只有6个（主对角线任意一侧的6个元素）。把B带入前面两个约束条件后可转化为：
上面两约束中的式子均可写为通式
的形式，定义3X3的单应矩阵H=[h1 h2 h3]的第i列列向量：
将如下表达式代入上述的约束单项式：
为了简化表达形式，令：
则有：
由此，两约束条件最终可以转化为如下形式：
如果我们拍摄了n张不同角度的标定图片，因为每张图片我们都可以得到一组（2个）上述的等式。其中，v12,v11,v22可以通过前面已经计算好的单应矩阵得到，因此是已知的，而b中6个元素是待求的未知数。因此，至少需要保证图片数 n>=3，我们才能解出b。
根据n张不同角度的标定图片，最终我们得到了一个矩阵集合 Vb=0 ，其中V是一个 (2n x 6) 的矩阵。如果 n>=3，就可以得到唯一解b（带有一个比例因子）。
如果 n=2，也就是只有两张标定图片，那么我们可以设置内参中的γ=0（γ表示由于制造误差产生的两个坐标轴偏斜参数，通常很小，可忽略），将前面式子（搬运到下图）中γ=0可以看到对应 B12=0，换句话说，就是增加了一个约束条件：[0, 1, 0, 0, 0, 0]b = 0。
如果n=1，只能假设u0, v0已知（位于图像中心）且 γ=0，这样只能解出fx, fy两个参数。
前面说到，B中包含一个尺度因子λ，即：
假设我们已经根据前面计算得到了矩阵B元素的值，那么根据已知的矩阵B很容易解出内参，如下：
得到内参数后，内参矩阵M也已知。单应矩阵H也已知，因此可继续求得外参数：
其中又由旋转矩阵性质有
则可得
实际情况下，数据中是存在噪音的，所以计算得到的旋转矩阵R并不一定能满足旋转矩阵的性质。所以通常根据奇异值分解来得到旋转矩阵R。
上述的推导结果是基于理想情况下的解，从理论上证明了张氏标定算法的可行性。但在实际标定过程中，一般使用最大似然估计进行优化。假设我们拍摄了n张标定图片，每张图片里有m个棋盘格角点。三维空间点X在图片上对应的二维像素为x，三维空间点经过相机内参M，外参R，t变换后得到的二维像素为x'，假设噪声是独立同分布的，我们通过最小化x, x'的位置来求解上述最大似然估计问题：
现在我们来考虑透镜畸变的影响，由于径向畸变的影响相对较明显，所以主要考虑径向畸变参数，根据经验，通常只考虑径向畸变的前两个参数k1,k2就可以（增加更多的参数会使得模型变的复杂且不稳定）。实际求解中，通常把k1,k2也作为参数加入上述函数一起进行优化，待优化函数如下所示
上述非线性优化问题通常用Levenberg-Marquardt（LM）算法进行迭代求解。一般将k1,k2初值设为0。
至此，张氏相机标定理论部分已经介绍完毕。