一、内参标定
相机内参主要标定相机的内参与畸变参数。
给标定板建立坐标系，建立归一化平面与标定板坐标系的转换关系。
归一化平面是指相机坐标系下 Z = 1 时的平面。
归一化平面上的点为：
(x, y, 1)
相机内参矩阵 K 为一个 3×3 矩阵：
fx 0 Cx
0 fy Cy
0 0 1
其中：
-
fx：X方向上的焦距，单位为像素 -
fy：Y方向上的焦距，单位为像素 -
Cx：相机主点在图像上的X坐标 -
Cy：相机主点在图像上的Y坐标
那么归一化平面坐标 (x, y) 与像素坐标 (u, v) 可以建立关系：
u = fx * x + Cx
v = fy * y + Cy
其中：
-
(x, y)是归一化平面坐标 -
(u, v)是图像上的像素坐标
设标定板坐标系下的某个点为：
(Xw, Yw, Zw, 1)
那么可以建立完整的投影关系：
s * (u, v, 1)T = K * [R, T] * (Xw, Yw, Zw, 1)T
其中：
-
K：相机内参 -
R：标定板坐标系到相机坐标系的旋转矩阵 -
T：标定板坐标系到相机坐标系的平移向量 -
s：齐次坐标的比例系数
将旋转矩阵写成：
R = (r1, r2, r3)
因为标定板是一个平面，所以可以直接把标定板坐标系建立在标定板平面上。
此时标定板上的所有角点都有：
Zw = 0
因此：
s * (u, v, 1)T
=
K * [r1, r2, r3, T] * (Xw, Yw, 0, 1)T
因为 Zw = 0，所以 r3 对应的这一项永远乘以0，因此可以删除。
得到：
s * (u, v, 1)T
=
K * [r1, r2, T] * (Xw, Yw, 1)T
定义：
H = K * [r1, r2, T]
于是：
s * (u, v, 1)T
=
H * (Xw, Yw, 1)T
这个 H 就叫做单应矩阵 Homography。
二、先求每一张标定图片的单应矩阵H
将 H 写成：
H =
h11 h12 h13
h21 h22 h23
h31 h32 h33
那么：
s * (u, v, 1)T
=
H * (Xw, Yw, 1)T
展开得到：
s*u = h11*Xw + h12*Yw + h13
s*v = h21*Xw + h22*Yw + h23
s = h31*Xw + h32*Yw + h33
把第三个公式：
s = h31*Xw + h32*Yw + h33
代入前面的两个公式。
可以得到：
u =
(h11*Xw + h12*Yw + h13)
/
(h31*Xw + h32*Yw + h33)
以及：
v =
(h21*Xw + h22*Yw + h23)
/
(h31*Xw + h32*Yw + h33)
将分母移到等式左边：
u*(h31*Xw + h32*Yw + h33)
=
h11*Xw + h12*Yw + h13
v*(h31*Xw + h32*Yw + h33)
=
h21*Xw + h22*Yw + h23
整理以后：
Xw*h11 + Yw*h12 + h13
-u*Xw*h31 - u*Yw*h32 - u*h33
= 0
Xw*h21 + Yw*h22 + h23
-v*Xw*h31 - v*Yw*h32 - v*h33
= 0
因此，一个标定板角点可以提供两个方程。
也就是说，只要我们知道：
标定板真实坐标 (Xw, Yw)
以及：
这个角点在图片中的像素坐标 (u, v)
就可以建立两个关于 H 的方程。
三、为什么可以求出H
H 一共有9个参数：
h11
h12
h13
h21
h22
h23
h31
h32
h33
但是单应矩阵是一个齐次矩阵。
例如：
H
和：
2H
所表达的投影关系是一样的。
因此 H 实际只有8个独立自由度。
一个角点能够提供两个方程。
所以理论上：
4个不共线的角点
就可以提供：
4 × 2 = 8
个约束，从而求出单应矩阵 H。
实际标定时，一个棋盘格一般会有几十个角点。
例如：
标定板坐标 图像像素坐标
(0, 0) -> (523.2, 314.5)
(20, 0) -> (557.8, 312.1)
(40, 0) -> (592.1, 309.8)
...
这些点全部参与计算。
由于实际角点检测存在误差，所以一般不会只使用4个点，而是使用所有角点，通过最小二乘或者SVD求一个最合理的 H。
第一张图片得到：
H1
第二张图片得到：
H2
第三张图片得到：
H3
最终得到：
H1
H2
H3
...
Hn
四、通过H计算相机内参K
前面已经知道：
H = K * [r1, r2, T]
将 H 按照列拆开：
H = (h1, h2, h3)
于是存在：
h1 = λ * K * r1
h2 = λ * K * r2
h3 = λ * K * T
其中 λ 是比例系数。
那么：
r1 = (1/λ) * K^-1 * h1
r2 = (1/λ) * K^-1 * h2
这里有一个非常关键的信息：
r1 和 r2 是旋转矩阵 R 的前两列。
而旋转矩阵的每一列其实都是一个单位坐标轴。
所以它们满足两个重要性质。
1. r1和r2互相垂直
因此：
r1T * r2 = 0
2. r1和r2长度相同
因为它们都是单位向量，所以：
r1T * r1 = r2T * r2 = 1
因此：
r1T * r1 = r2T * r2
这两个条件就是后面求相机内参的关键。
五、利用旋转矩阵的性质求K
因为：
r1 = (1/λ) * K^-1 * h1
r2 = (1/λ) * K^-1 * h2
根据：
r1T * r2 = 0
代入得到：
h1T * K^-T * K^-1 * h2 = 0
再根据：
r1T * r1 = r2T * r2
可以得到：
h1T * K^-T * K^-1 * h1
=
h2T * K^-T * K^-1 * h2
为了让公式简单一点，定义：
B = K^-T * K^-1
那么上面的两个公式就变成：
h1T * B * h2 = 0
以及：
h1T * B * h1
-
h2T * B * h2
=
0
所以：
每一张标定板图片，都可以利用自己的单应矩阵
H提供两个关于相机内参的约束。
六、B矩阵是什么
因为：
B = K^-T * K^-1
所以 B 是一个对称矩阵。
可以写成：
B11 B12 B13
B12 B22 B23
B13 B23 B33
因此它只有6个不同的未知数：
B11
B12
B22
B13
B23
B33
把它们写成一个向量：
b =
(B11,
B12,
B22,
B13,
B23,
B33)T
每张标定图片可以提供两个约束。
假设拍摄了很多张不同角度的标定图片，就可以得到很多方程。
最后组合成：
V * b = 0
然后通过SVD计算出：
b
从而得到矩阵：
B
七、通过B恢复内参K
因为已经知道：
B = K^-T * K^-1
而内参矩阵为：
fx 0 Cx
0 fy Cy
0 0 1
所以可以通过 B 反推出：
fx
fy
Cx
Cy
最终得到：
K =
fx 0 Cx
0 fy Cy
0 0 1
如果考虑像素坐标轴之间不完全垂直，还存在一个 skew 参数：
K =
fx skew Cx
0 fy Cy
0 0 1
但是现代工业相机通常认为：
skew = 0
所以平时看到的基本都是：
fx 0 Cx
0 fy Cy
0 0 1
到这里，已经得到了相机内参的初始值。
八、为什么还要计算每张图片的R和T
前面：
H = K * [r1, r2, T]
现在 K 已经知道了，所以可以反过来计算每张标定图片对应的外参。
根据：
h1 = λ * K * r1
可以得到：
r1 = λ * K^-1 * h1
同理：
r2 = λ * K^-1 * h2
以及：
T = λ * K^-1 * h3
比例系数 λ 可以通过下面的关系计算：
λ = 1 / ||K^-1 * h1||
求出 r1 和 r2 后：
r3 = r1 × r2
于是得到：
R = [r1, r2, r3]
最终每一张图片都可以得到一套：
R
T
需要注意：
同一个相机的内参K是不变的。
但是每一张标定图片的：
R
T
是不一样的。
因为拍照的时候，标定板相对于相机的位置和姿态一直在变化。
例如：
第1张图片：
R1
T1
第2张图片：
R2
T2
第3张图片：
R3
T3
而所有图片共享同一个：
K
九、为什么还要计算畸变
前面使用的模型：
u = fx*x + Cx
v = fy*y + Cy
是一个理想针孔相机模型。
但是实际相机里面有镜头。
真实镜头并不是完全理想的，所以图像会发生畸变。
例如原本应该是一条直线：
-------------------------
可能拍出来变成：
----_____________----
或者：
____-------------____
因此只计算：
fx
fy
Cx
Cy
还不够。
还需要计算镜头的畸变参数。
常见畸变主要分为：
- 径向畸变
- 切向畸变
十、径向畸变
假设理想归一化平面坐标为：
(x, y)
先计算这个点距离归一化平面中心的距离：
r² = x² + y²
考虑径向畸变以后：
xd =
x * (1 + k1*r² + k2*r⁴ + k3*r⁶)
yd =
y * (1 + k1*r² + k2*r⁴ + k3*r⁶)
其中：
k1
k2
k3
就是径向畸变参数。
径向畸变最大的特点就是：
距离图像中心越远，畸变通常越明显。
常见的桶形畸变、枕形畸变，主要就是径向畸变。
十一、切向畸变
如果镜头和相机内部的成像平面没有做到完全平行，还会产生切向畸变。
切向畸变通常使用：
p1
p2
描述。
公式为：
xd =
x
+ 2*p1*x*y
+ p2*(r² + 2*x²)
yd =
y
+ p1*(r² + 2*y²)
+ 2*p2*x*y
十二、径向畸变和切向畸变合在一起
实际计算的时候，通常会把两种畸变一起计算。
得到：
xd =
x*(1 + k1*r² + k2*r⁴ + k3*r⁶)
+ 2*p1*x*y
+ p2*(r² + 2*x²)
yd =
y*(1 + k1*r² + k2*r⁴ + k3*r⁶)
+ p1*(r² + 2*y²)
+ 2*p2*x*y
得到畸变后的归一化坐标：
(xd, yd)
再通过相机内参转成像素：
u = fx*xd + Cx
v = fy*yd + Cy
所以OpenCV中常见的5个畸变参数就是：
k1
k2
p1
p2
k3
十三、畸变参数到底是怎么算出来的
现在我们已经知道：
- 标定板角点真实坐标
- 每个角点真实拍摄到的像素坐标
- 初步计算出来的相机内参K
- 每张图片对应的R和T
假设标定板中一个角点为：
Pw = (Xw, Yw, Zw)
首先通过当前图片对应的外参：
Pc = R * Pw + T
得到这个点在相机坐标系中的坐标：
Pc = (Xc, Yc, Zc)
然后投影到归一化平面：
x = Xc / Zc
y = Yc / Zc
然后加入当前的畸变参数：
(x, y)
↓
畸变模型
↓
(xd, yd)
然后通过内参得到预测像素：
u' = fx*xd + Cx
v' = fy*yd + Cy
也就是说，根据当前相机参数，我们预测这个角点应该出现在：
(u', v')
但是实际图片检测出来的位置可能是：
(u, v)
那么二者之间就存在误差：
du = u - u'
dv = v - v'
这个点的误差大小为：
e = sqrt(du² + dv²)
这个就是重投影误差。
十四、为什么叫重投影误差
因为整个过程实际上是：
已知标定板3D坐标
↓
R、T
↓
相机坐标系3D坐标
↓
/ Z
↓
归一化平面坐标
↓
畸变模型
↓
K
↓
预测像素坐标
然后拿这个预测出来的像素坐标，和真实图片里面检测到的像素坐标进行比较。
相当于：
把标定板上的三维点重新投影回图片上。
所以叫：
Reprojection Error，重投影误差。
十五、最终不是只计算一次，而是整体优化
前面通过单应矩阵计算出来的：
fx
fy
Cx
Cy
以及：
R
T
主要是为了得到一个比较好的初始值。
最终标定还需要做非线性优化。
优化变量包括：
fx
fy
Cx
Cy
k1
k2
k3
p1
p2
第1张图片的 R1、T1
第2张图片的 R2、T2
第3张图片的 R3、T3
...
然后计算所有图片、所有角点的重投影误差。
目标就是：
让所有角点的预测像素
尽可能接近
真实检测出来的像素
也就是最小化：
Σ Σ || 实际像素 - 预测像素 ||²
不断调整：
K
畸变参数
每张图片的R、T
让整体误差越来越小。
最后得到最优的：
fx
fy
Cx
Cy
k1
k2
p1
p2
k3
这就是最终保存下来的相机标定参数。
十六、内参标定完整流程
整个内参标定过程可以总结成：
已知标定板尺寸
↓
建立标定板坐标系
↓
得到角点真实坐标 (Xw, Yw, 0)
↓
拍摄多张不同姿态的标定板图片
↓
检测图片角点 (u, v)
↓
建立：
(Xw, Yw) <-> (u, v)
对应关系
↓
每张图片计算单应矩阵 H
↓
利用旋转矩阵正交约束
计算相机内参 K
↓
计算每张图片对应的 R、T
↓
加入镜头畸变模型
↓
计算重投影误差
↓
对 K、畸变、R、T 进行整体非线性优化
↓
得到最终标定结果
最终得到：
相机内参：
fx
fy
Cx
Cy
以及：
畸变参数：
k1
k2
p1
p2
k3
十七、一句话理解内参标定
内参标定本质上就是：
我知道标定板上的一个点在真实世界中的位置，也知道这个点最终拍到了图片的哪个像素，通过大量这样的对应关系，反推出这个相机到底是如何完成“3D世界 -> 2D图像”这个投影过程的。
最终得到两类参数：
K：描述理想针孔相机模型
D：描述真实镜头产生的畸变
也就是：
真实3D点
↓
外参 R、T
↓
相机坐标
↓
透视投影
↓
归一化坐标
↓
畸变
↓
相机内参 K
↓
像素坐标
然后通过不断减小预测像素和实际像素之间的重投影误差，最终得到最符合真实相机的内参与畸变参数。
