【Python计算机视觉】相机标定(Camera Calibration)

一、相机投影

1.1 针孔相机模型

1.用一个带有小孔的板遮挡在屏幕与物之间，屏幕上就会形成物的倒像，我们把这样的现象叫小孔成像。前后移动中间的板，屏幕上像的大小也会随之发生变化，这种现象说明了光沿直线传播的性质。

2.相机将三维世界中的坐标点（单位：米）映射到二维图像平面（单位：像素）的过程能够用一个几何模型来描述，其中最简单的称为针孔相机模型 (pinhole camera model)，其框架如下图所示。

很显然，坐标的转换便是一个降维的过程。为了确定这种三维到二维的映射关系，我们进行如下对原理的讨论。

1.2 原理

1.2.1 坐标系标定

1.世界坐标系（world coordinate）(xw,yw,zw) ，也称为测量坐标系，是一个三维直角坐标系，以其为基准可以描述相机和待测物体的空间位置。世界坐标系的位置可以根据实际情况自由确定。
2.相机坐标系（camera coordinate）(xc,yc,zc)，也是一个三维直角坐标系，原点位于镜头光心处，x、y轴分别与相面的两边平行，z轴为镜头光轴，与像平面垂直。
3.像素坐标系 uov是一个二维直角坐标系，反映了相机中像素的排列情况。原点o位于图像的左上角，u轴、v轴分别于像面的两边平行。像素坐标系中坐标轴的单位是像素（整数）。
4.图像坐标系 XOY，其坐标轴的单位通常为毫米（mm），原点是相机光轴与相面的交点（称为主点），即图像的中心点，X轴、Y轴分别与u轴、v轴平行。故两个坐标系实际是平移关系，即可以通过平移就可得到。
5.世界坐标系到相机坐标系的标定：
第一步是从世界坐标系到相机坐标系的变换，这一步是三维点到三维点的变换，属于刚体变换：即物体不会发生形变，只需要进行旋转和平移。
R：表示旋转矩阵（相机外参）
T：表示偏移向量（相机外参）

6.相机坐标系到图像坐标系的标定：
第二步是相机坐标系转为图像坐标系，这一步是三维点到二维点的变换，即从3D转换到2D.属于透视投影关系。

得到矩阵公式：

由于此时投影点p的单位还是mm，并不是像素pixel，需要进一步转换到像素坐标系。

7.图像坐标系到像素坐标系的标定：
像素坐标系和图像坐标系都在成像平面上，只是各自的原点和度量单位不一样。图像坐标系的原点为相机光轴与成像平面的交点，通常情况下是成像平面的中点或者叫principal point。图像坐标系的单位为mm，属于物理单位，而像素坐标系的单位是pixel，我们平常描述一个像素点都是几行几列。所以这两者之间的转换如下：
其中dx和dy表示每一列和每一行分别代表多少mm，即1pixel = dx mm

于是由以上四个坐标系的转换得到如下式子：

右式结果的第一个矩阵为相机内参，第二个矩阵为相机外参，可由张有正标定法获取（戳这里查看张有正标定有关信息）

1.2.2 相机畸变现象影响

由于光学透镜固有的透视失真导致相机镜头畸变的现象我们统称为相机畸变。此处我们着重讨论图像径向畸变：
产生原因：

• 透镜质量原因
• 光线在远离透镜中心的地方比靠近中心的地方更加弯曲

畸变的本质是光学系统在不同视场由于镜头屈光度与光阑位置不同引起像高放大率不同导致的。

分类：
1.桶状畸变：又称“桶形失真”，视野中光轴中心区域的放大率远大于边缘区域的放大率，常出现在广角镜头和鱼眼镜头中。
2.枕形畸变：又称“鞍形形变”，视野中边缘区域的放大率远大于光轴中心区域的放大率，常用在远摄镜头中。

矫正：

调节公式为：

• 公式里（x0，y0）是畸变点在成像仪上的原始位置，（x，y）是校正后新的位置。
• 实现过程是，对输出图的点做遍历——以1080p的图像为例，从点（0,0）到点（1919,1079），一行一行的遍历——依次找到输出点（x,y）对应的原图点（x0, y0）的像素值，再将（x0, y0）的值赋给（x, y）。如果计算出来的对应的原图的点（x0,y0）不是整数，则用二次线性插值计算此点，然后赋值给（x, y）。

1.2.3 相机旋转平移影响

由于日常生活中在进行拍摄时，相机的位置是不断变化的，所以会有相机旋转与平移对拍摄的影响。所以将上面坐标转换得到的x ~ K[I|0]X矩阵进行修改，得到x ~ K[R|t]X矩阵（外部参数矩阵）：

1.相机旋转包括三个方向：X,Y,Z，（即pitch，yaw，roll）

2.将三个方向上的旋转相乘，得到旋转矩阵R：

1.2.4 两类重要参数

二、相机标定

同步标定内部参数和外部参数，一般包括两种策略s：

1. 光学标定：利用已知的几何信息（如定长棋盘格）实现参数求解。
2. 自标定：在静态场景中利用structure from motion估算参数。

2.1 线性回归（最小二乘）

线性回归（最小二乘）单单考虑坐标的转换：

通过空间中已知的坐标的（特征）点（Xi，Yi，Zi），以及它们在图像中的对应坐标（ui，vi），直接估算11个待求解的内部和外部参数。

表示为矩阵形式：

即可用最小二乘法求解。
（若要回顾单硬性变换，请点击这里）

给定超定方程：Ax = 0 后，x的解即为A^TA最小特征值对应的特征向量。

标定参数线性回归的优缺点：
优点：

• 所有相机参数集中在一个矩阵中，便于求解
• 所有矩阵可以直接描述世界坐标中的三维点，到二维图像平面中点的映射关系。

缺点：

• 无法直接得知具体的内参数和外参数
• 求解出的11个未知量，比待标定参数（9个）更多。带来了参数不独立/相关问题
• 对噪声/误差敏感。

2.2 非线性优化

用概率的角度去看最小二乘问题。

• 特征点的投影方程为：
  
• 给定{（ui，vi）}，标定参数矩阵M的概率为：
  
• 给定{（ui，vi）}，标定参数矩阵M的似然函数为：
  
• 相应求解策略：牛顿方法、高斯-牛顿方法、Levenberg-Marquardt算法等。

2.3 张有正标定

（1）计算单应性矩阵H

设三维世界坐标的点为X=[X,Y,Z,1]^T
二维相机平面像素坐标为m=[u,v,1]^T
所以标定用的棋盘格平面得到图像平面的单应性关系为：
其中，s为尺度因子，K为摄像机内参数，R是旋转矩阵，t为平移向量，令：

s做为齐次坐标来说，不会改变齐次坐标值。将世界坐标系构造在棋盘格平面上，令棋盘格平面为Z=0的平面。则：
将K[r1,r2,t]叫做单应性矩阵H，即：

H是一个齐次矩阵，所以有八个未知数，至少需要八个方程，每对对应点能提供两个方程，所以至少需要四个对应点，就可以算出世界平面到图像平面的单应性矩阵H。

（2）计算内参数矩阵


可以看出，B是一个对称阵，所以B的有效元素为6个，让这6个元素写成向量b，即：


利用约束条件可以得到：

通过上式，我们至少需要三幅包含棋盘格的图像，可以计算得到B，然后通过cholesky分解，得到相机的内参数矩阵K。
（3）计算外参数矩阵
由之前的推导可得：

（4）最大似然估计
由于可能存在高斯噪声，所以使用最大似然估计进行优化。这里使用的是多参数非线性系统优化问题的Levenberg-Marquardt算法进行迭代优化求最优解。

（5）径向畸变估计
数学表达式：

化作矩阵形式：

得到：

由此计算得到畸变系数k
综上，我们得到了相机的内参、外参和畸变系数。

2.4 相机标定实现

2.4.1 算法描述

1. 打印一张棋盘格A4纸张（黑白间距已知），并贴在一个平板上