Opencv2.4.9源码分析——Stitching（三）

3、相机参数评估

3.1 原理

相机参数的评估也称为相机定标。要想理解这部分内容，首先应该从成像原理开始讲起。

图6 小孔成像原理

从图6可以看出，真实物体通过小孔映射到成像平面上，小孔到成像平面的距离称为焦距f。在成像平面上的图像是镜像倒立的，所以为了研究方便，在小孔和物体之间定义一个虚拟成像平面（在后面，我们把该平面也称为成像平面），它与小孔的距离也为焦距，则两个成像平面的图像大小是相同的，但虚拟成像平面上的图像与原物体的方向是一样的。

图7 成像的几何模型

我们以小孔为坐标原点建立一个三维直角坐标系XYZ（如图7所示），坐标原点C称为相机的光心。成像平面xy平行于XY，并且距离光心C为f，其中Z轴定义为光轴，它与成像平面xy的交点为P，因此CP＝f。设空间中的任一点Q的坐标为（X, Y, Z），该点映射到成像平面的点q的坐标为(x, y,f)。

由几何知识可得：

（31）

式中，λ＝Z/f，则x＝fX/Z，y＝fY/Z。因此空间中Q点的三维坐标映射到成像平面的二维坐标q点在齐次坐标下的线性映射关系为：

（32）

如图7所示，在成像平面（坐标系为xy）的坐标原点为P，但图像（设坐标系为uv）的坐标原点一般在左上角，所以这两个坐标系之间需要通过平移来进行转换。另外成像平面xy是长度单位，而相机图像传感器的单位是像素，因此像素与长度之间也是需要转换的，而且水平和垂直的像素往往是不相同的，所以横纵轴的转换系数是不一致。因此式32改写为：

（33）

式中，f_u和f_v为焦距f在横纵轴的长度和像素的转换，它们之间的关系可以写为f_v＝αf_u，(c_x, c_y)为坐标平移。因为矩阵K的参数都是基于相机内部自身的参数，因此K称为相机内参数，K＝TV。

相机除了具有内参数外，还有外参数。式33中的(X, Y, Z)称为相机坐标，而它与真实的世界坐标(X’, Y’, Z’)还存在欧几里得变换，即：

（34）

式中，R为3×3的旋转矩阵，t为3×1的平移向量。式34代入式33，得

（35）

式中，[R|t]称为相机外参数，M称为投影矩阵。

如果得到的两幅图像如图3所示的那样，即相机在三角架上通过旋转得到的两幅图像，则对于同一个世界坐标上的点(X’, Y’, Z’)，在两幅图像上的坐标点分别为(u₀, v₀)和(u₁, v₁)，即：

（36）

由于相机只做旋转处理（或者我们认为物体离相机很远），则t₀＝t₁＝0，从而得到(u₀, v₀)和(u₁, v₁)的关系为：

（37）

式中，R₁₀为由图像1到图像0的相对旋转矩阵，R₁₀＝R₁R₀^-1。由式2可得：

（38）

为简单起见，我们设两幅图像的相机内参数中的坐标平移都为0，即K＝V，因为抛弃参数T对图像拼接影响不大。

我们在评估焦距时，还需定义f_1u＝f_1v＝f₁，f_0u＝f_0v＝f₀，即设图像的长宽等像素：

（39）

则式38表示为

（40）

式中，R₁₀＝[r_ij]。

由式40我们就可以得到焦距f₀和f₁：观察矩阵R₁₀可知，R₁₀前两行一定有相同的范数，并且是正交的，因此

（41）

（42）

由式41可得：

（43）

由式42可得：

（44）

由式43和式44得到了两个f₀，选取哪个呢？比较式43和式44中分母部分的绝对值的大小，如果式43的分母大，则选择分式大的值作为f₀，否则如果式44的分母大，则选择值小的作为f₀。

同理，矩阵R₁₀的前两列也一定有相同的范数，并且也是正交的，因此

（45）

（46）

由式45可得

（47）

由式46可得

（48）

比较式47和式48中分母部分的绝对值的大小，如果式48的分母大，则选择分式大的值作为f₁，否则如果式47的分母大，则选择值小的作为f₁。

如果两幅图像的焦距相同，则最终这两幅图像的焦距f为：

（49）

当评估计算得到多个f时，可取这些f的中值作为所有相机的焦距。

焦距得到后，我们就可以由式38得到R₁₀：

（50）

则R₀₁（由图像0到图像1的相对旋转矩阵）为：

（51）

则

（52）

对于刚性物体，它的旋转都是沿笛卡尔坐标系的x轴、y轴和z轴旋转，则分别沿着这三个轴的旋转矩阵定义为：

（53）

则旋转矩阵R表示为：

（54）

三维旋转除了可以用旋转矩阵描述外，还可以用旋转向量r描述，即r＝[r_x, r_y, r_z]^T。旋转向量的长度（模）表示绕轴逆时针旋转的角度θ。旋转向量和旋转矩阵可以通过Rodrigues算法进行转换。由旋转向量转换为旋转矩阵的Rodrigues算法描述如下：

（55）

（56）

（57）

式中，I为3×3的单位矩阵。而由旋转矩阵转换为旋转向量的Rodrigues算法公式为：

（58）

当有多幅图像需要拼接为一幅图像时，是要以其中一幅图像为基准，其他图像都要旋转到该基准图像平面上的，所以就需要找到基准平面。这里用到的算法为最大生成树算法。

待拼接图像的排列是无序的，而且我们事先是不知道它们之间的关系的，我们只知道它们之间的单应矩阵，而单应矩阵是由图像间的内点计算评估得到的。由此我们可以构造一个无向图G，G的节点为图像，G的边为内点数，然后利用并查集在该G中得到一棵最大生成树。

图8 最大生成树

图8为用于拼接的最大生成树的一个例子，图(a)为无向图，节点为图像（A、B、C、D、E），节点间的边为内点数。图(b)为最大生成树，由图像C到图像B要经过最大的边连接，所以要经过图像A，而图像C和图像B之间的连接就需要去掉了。

我们把树的中心节点作为基准图像。中心节点的确定方法为：计算每一个节点到所有叶节点的距离，把其中的最大值作为该节点的值；然后选择这些值中最小者作为中心节点。这里的距离指的是节点间的节点数。如图8(c)所示，节点A和C为中心节点。中心节点可能是1个，也可能是2个，如果是2个，则选择其中一个即可。

基于以上方法，我们得到了相机的内外参数，但这样得到的参数忽略了多个图像间的约束，而且会产生累计误差。这时，我们就需要用到光束平差法（Bundle Adjustment）来精确化相机参数。光束指的是相机光心“发射”出来的光束（或射线），它透过相片达到物点，因此相片中的点应该和物点处于一个光束线上，但当两者不共线时，我们就需要对结构和视角参数进行调节，以达到最优解甚至共线的目的。最优化一般采用前文介绍过的LM算法。

应用于光束平差法的LM算法，误差指标函数可以有两个，一个是重映射误差，另一个是射线发散误差。

重映射误差公式为式25（即一个内点要有x轴和y轴两个误差值），而单应矩阵H是由式38得到。也就是说H是由相机的内外参数得到。相机的内外参数一共有7个：f_u、α、c_x、c_y、r_x、r_y和r_z。前4个参数是内参数（见式33），后3个参数是外参数（即式55中的旋转向量的三个元素）。因此式25中的h为h(f_u, α, c_x, c_y, r_x, r_y,r_z)，由此得到J(h)为：

（59）

式中，n表示待拼接的图像数量，也就是所有的相机。有时为了计算方便，导数可以用差分近似，例如我们要计算e对f_u的偏导，则

（60）

式中，∆表示一个很小的数。

图9 射线发散概念

第二种误差指标函数是基于射线发散原理。如图9所示，不同的相机发出的射线透过相片后达到同一个物点，但由于误差，两条射线不会相交，或者称为两条射线发散了，我们就把这两条射线间的最短距离d定义为射线发散误差。

下面，我们不加证明的给出射线发散误差的计算公式：

设(u₁, v₁, 1)和(u₂, v₂, 1)为两幅不同图像的同一特征点的齐次坐标，则由单应矩阵H₁和H₂分别得到它们的物点坐标为：

（61）

分别对坐标进行归一化处理：

（62）

（63）

为了简化计算，射线间的最短距离可以表示为基于不同焦距f₁和f₂下的归一化后的物点坐标之差：

（64）

式64说明在计算射线发散误差时，每一个内点有x轴、y轴和z轴三个误差值。

我们假设射线就是相机的光轴，因此单应矩阵H中的参数α＝1，c_x＝c_y＝0，所以式25中h只是基于4个参数的变量，即h(f_u, r_x,r_y, r_z)，由此得到射线发散误差的J(h)为：

（65）

前面介绍的光束平差法会引起波形效应，即拼接的图像会呈现蛇形分布，这是因为真实拍摄相片时不大可能都保持水平而不倾斜的，也就是重力轴没有垂直于图像。因此我们就需要引入一个全局校正矩阵Q，用于“拉直”拼接图像，该方法也成为波形校正。

校正的目的应该是，在第k个相机旋转矩阵R_k乘以Q后，全局y轴应该与图像的x轴垂直，这种约束条件可以表示为：

（66）

式中，i＝(1, 0, 0)，j＝(0, 1, 0)，该式的含义是R_k的第一行r_k0垂直于Q的第二列q₁。式66的这类约束问题可以看成是最小二乘问题：

（67）

因此，q₁是矩量矩阵∑r^T