【相机标定与三维重建原理及实现】学习笔记2——标定的算法与方法
作者：互联网
目录
前言
由于本人近期正在展开数字图像相关技术用于测量材料形变方向的研究，既然涉及到使用图像处理参与到测量或检测研究当中，就肯定避不开构建物空间上的任意一点与相机所拍摄到的图像上一点之间的数学关系及数学模型。之前有过标定经验的我只是单纯的使用别人封装好的函数，如今仔细推导过后才发现其中奥秘所在。我认为想要在一个技术上有所创新，最核心的还是能够将其学的透彻，我希望将自己作为一个初学者学习和推导该原理的过程记录下来，也方便之后每一个涉足该领域知识的人能更清晰更快的应用这些知识。
本文所写的内容主要参考《学习OpenCV 3 》1以及 大奥特曼打小怪兽大佬的博客第六节、双目视觉之相机标定2。但是本文的讲解思路会稍有不同，同时也会更偏向数学推导。如果本文对你的帮助不大，建议可以看一下提到的书籍与博客，可以提供一些不同角度的讲解。
本文的推导基于上一篇笔记【相机标定与三维重建原理及实现】学习笔记1——相机模型数学推导详解，建议先行阅读，上一篇笔记探讨了为什么要标定以及标定是在标定什么，本章节我希望探讨的有以下几个问题：
- 如何进行标定？
- 标定时候究竟需要多少标定图像？每个标定图像需要有多少特征点数？
本专栏目录：
在上一篇笔记中，我们知道了标定的主要目的是为了建立摄像机图像像素位置与物体空间位置之间的关系，即世界坐标系与像素坐标系之间的关系。而这些关系里面主要有两种模型需要解算：小孔成像模型与畸变模型，其中前者是线性的，后者是非线性的。首先来看一下小孔成像模型的表达的映射关系：
z
c
[
u
v
1
]
=
[
f
x
0
c
x
0
f
y
c
y
0
0
1
]
[
R
t
0
3
T
1
]
[
x
w
x
w
z
w
1
]
z_{c}\left[\begin{array}{l} u \\ v \\ 1 \end{array}\right] =\left[\begin{array}{ccc} f_{x} & 0 & c_{x}\\ 0 &f_{y}& c_{y}\\ 0 & 0 & 1 \end{array}\right]\left[\begin{array}{cccc} R & t \\ 0_{3}^{T} & 1 \end{array}\right]\left[\begin{array}{c} x_{w} \\ x_{w} \\ z_{w} \\ 1 \end{array}\right]
zc⎣⎡uv1⎦⎤=⎣⎡fx000fy0cxcy1⎦⎤[R03Tt1]⎣⎢⎢⎡xwxwzw1⎦⎥⎥⎤其中
- u , v u,v u,v表示像素坐标系下的坐标
- f x = f ⋅ s x ， f y = f ⋅ s y f_{x}=f \cdot s_{x}，f_{y}=f \cdot s_{y} fx=f⋅sx，fy=f⋅sy为两个组合量，表达摄像机系统的焦距与对应方向的单位尺寸像素数相乘的结果
- R R R表示从世界坐标系到相机坐标系的旋转矩阵，其实际包含的未知量只有三个，即坐标系绕三个轴的旋转角度
- t t t表示从世界坐标系到相机坐标系的平移向量，包含三个未知量
- x w , x w , z w x_{w},x_{w},z_{w} xw,xw,zw表示对应物点在世界坐标系下的坐标位置
- z c z_{c} zc表示物点在相机坐标系下的Z方向坐标值，其具体数值可以从 [ R t 0 3 T 1 ] [ x w x w z w 1 ] \left[\begin{array}{cccc} R & t \\ 0_{3}^{T} & 1 \end{array}\right]\left[\begin{array}{c} x_{w} \\ x_{w} \\ z_{w} \\ 1 \end{array}\right] [R03Tt1]⎣⎢⎢⎡xwxwzw1⎦⎥⎥⎤中提取出来，是一个中间变量
在计算机视觉中将这种从一个平面到另一个平面的投影映射称之为单应性，我们将上面这个公式稍微改写一下，将
z
c
z_{c}
zc移动到等式的右边，设其倒数为
s
=
1
/
z
c
s=1/z_{c}
s=1/zc，设摄像机系统的内参矩阵为
M
=
[
f
x
0
c
x
0
f
y
c
y
0
0
1
]
M=\left[\begin{array}{ccc} f_{x} & 0 & c_{x}\\ 0 &f_{y}& c_{y}\\ 0 & 0 & 1 \end{array}\right]
M=⎣⎡fx000fy0cxcy1⎦⎤，外参矩阵为
W
=
[
R
t
0
3
T
1
]
W=\left[\begin{array}{cccc} R & t \\ 0_{3}^{T} & 1 \end{array}\right]
W=[R03Tt1]，则我们可以得到变化后的公式：
[
u
v
1
]
=
s
⋅
M
⋅
W
[
x
w
x
w
z
w
1
]
\left[\begin{array}{l} u \\ v \\ 1 \end{array}\right]=s\cdot M\cdot W\left[\begin{array}{c} x_{w} \\ x_{w} \\ z_{w} \\ 1 \end{array}\right]
⎣⎡uv1⎦⎤=s⋅M⋅W⎣⎢⎢⎡xwxwzw1⎦⎥⎥⎤设矩阵
H
=
M
⋅
W
H=M\cdot W
H=M⋅W，则有:
[
u
v
1
]
=
s
⋅
H
[
x
w
x
w
z
w
1
]
\left[\begin{array}{l} u \\ v \\ 1 \end{array}\right]=s\cdot H\left[\begin{array}{c} x_{w} \\ x_{w} \\ z_{w} \\ 1 \end{array}\right]
⎣⎡uv1⎦⎤=s⋅H⎣⎢⎢⎡xwxwzw1⎦⎥⎥⎤其中
H
H
H又被称为单应性矩阵，表达了两个坐标系之间的映射关系，而
s
s
s是H里面可以分解得到的一个变量（上面提到了
z
c
z_{c}
zc是一个中间变量，因此是可以从
H
H
H中提取的）。写成这样的形式，我想很容易就可以看出，这不就是个形如线性方程组的形式嘛！！那按照求解线性方程组的方式，只要知道足够多的相互对应的物点和像点，就可以构建大于等于未知量数目的方程组，这不就可以求解未知量了嘛！
实际上，标定正是通过这样的方式来求取图像与世界之间的映射关系的，但这个式子其实是一个非线性方程组（因为有
s
s
s的存在），如果要求解映射关系里的未知量就需要用到非线性最小二乘法，来找到一组最优的解使得目标函数（误差的平方和）最小，在这个问题中目标函数可以写为：
∑
i
=
1
n
∑
j
=
1
m
∥
m
i
j
−
s
⋅
H
⋅
w
i
j
∥
2
\sum_{i=1}^{n}\sum_{j=1}^{m}\left \| \mathbf{m}_{ij}- s\cdot H\cdot\mathbf{w}_{ij}\right \|^{2}
i=1∑nj=1∑m∥mij−s⋅H⋅wij∥2其中
m
i
j
\mathbf{m}_{ij}
mij表征了上述式子中的像点
[
u
v
1
]
\left[\begin{array}{l}u \\v \\1\end{array}\right]
⎣⎡uv1⎦⎤，
w
i
j
\mathbf{w}_{ij}
wij表征了上述式子中的物点
[
x
w
x
w
z
w
1
]
\left[\begin{array}{c} x_{w} \\ x_{w} \\ z_{w} \\ 1 \end{array}\right]
⎣⎢⎢⎡xwxwzw1⎦⎥⎥⎤，整个式子的物理意义即计算三维空间中的点
w
i
j
\mathbf{w}_{ij}
wij经过坐标映射投影回图像后的坐标与实际图像上提取的点
m
i
j
\mathbf{m}_{ij}
mij的坐标之间的2-范数（即两点间的距离），数值越小则表明误差越小、效果越好。 以这个目标函数为准则，利用非线性最小二乘即可求解出最优的单应性矩阵
H
H
H。而这正是标定的算法核心，但照这样直接去做标定存在两个问题：
- 被拍摄的物点在世界坐标系下的坐标怎么进行度量？
- 目标函数所包含的未知量很多，在初始值不好的情况（例如像点特征提取效果差或者物点的度量误差大）下，很容易陷入局部最优中。
而张正友大佬（膜拜！）提出的使用标定板进行标定的方法，很大程度上解决了这两个问题。
一.张氏标定法
“张氏标定”是指张正友教授于1998年提出的单平面棋盘格的摄像机标定方法。张氏标定法已经作为工具箱或封装好的函数被广泛应用。张氏标定的原文为“A Flexible New Technique forCamera Calibration”3。张正友大佬的这个方法其实听起来会觉得非常简单，他的方法有两个核心，第一，使用如下图所示的棋盘格标定板方便提取像点和衡量物点；第二，默认标定板平面为世界坐标系中
z
=
0
z=0
z=0的平面。 然后通过拍拍摄若干张这样的图片，计算出不同标定板姿态下的单应性矩阵，再利用多组单应性矩阵构成的超定方程组求解内参和外参。
这样一来，标志物方便了从图像上获取特征点；其次将标定板认为是
z
=
0
z=0
z=0的平面，只需要知道特征点打印的间距即可很方便地度量特征点在世界坐标系下的坐标位置（算法中通常会把四个角上提取到的特征点之一作为世界坐标系的原点）；最后由于所有物点都位于
z
=
0
z=0
z=0的平面上，这样一来大大减少了目标函数的待优化量，之前的映射公式可以重新改为：
[
u
v
1
]
=
s
⋅
M
⋅
W
[
x
w
x
w
z
w
1
]
=
s
⋅
M
⋅
[
r
1
r
2
r
3
t
]
[
x
w
x
w
0
1
]
=
s
⋅
M
⋅
[
r
1
r
2
t
]
[
x
w
x
w
1
]
\begin{aligned}\left[\begin{array}{l} u \\ v \\ 1 \end{array}\right]&=s\cdot M\cdot W\left[\begin{array}{c} x_{w} \\ x_{w} \\ z_{w} \\ 1 \end{array}\right]\\&=s\cdot M\cdot\begin{bmatrix} \mathbf{r}_{1} & \mathbf{r}_{2} & \mathbf{r}_{3} & \mathbf{t} \end{bmatrix} \left[\begin{array}{c} x_{w} \\ x_{w} \\ 0 \\ 1 \end{array}\right] \\ &=s\cdot M\cdot\begin{bmatrix} \mathbf{r}_{1} & \mathbf{r}_{2} & \mathbf{t} \end{bmatrix} \left[\begin{array}{c} x_{w} \\ x_{w} \\ 1 \end{array}\right] \end{aligned}
⎣⎡uv1⎦⎤=s⋅M⋅W⎣⎢⎢⎡xwxwzw1⎦⎥⎥⎤=s⋅M⋅[r1r2r3t]⎣⎢⎢⎡xwxw01⎦⎥⎥⎤=s⋅M⋅[r1r2t]⎣⎡xwxw1⎦⎤此时的
H
=
M
⋅
[
r
1
r
2
t
]
H=M\cdot\begin{bmatrix} \mathbf{r}_{1} & \mathbf{r}_{2} & \mathbf{t} \end{bmatrix}
H=M⋅[r1r2t]为一个3X3的矩阵，计算量相比于之前减少了很多。而正是这样的设计大大减少了算法的耗时并且提升了标定的鲁棒性！接下来我们需要探讨使用张氏标定法究竟需要拍多少张这样的图片以及标定板上究竟需要多少个特征点。
二.张氏标定需要的特征点数以及拍摄图片数
对张氏标定法的算法思路进行总结，可以得到如图所示的流程示意，使用单张标定图像时，在特征点足够的情况下，我们可可以计算得到在当前标定板姿态（外参）的情况下的单应性矩阵 H H H，但这并不足以我们分解得到摄像机系统的内参，因为 H H H是外参和内参的组合量，如果要计算其中内参部分的具体数值我们还需要更多的方程，因此我们需要拍摄不同标定板姿态下的图像。换而言之，对于张氏标定法具体需要的标定板特征点数和拍摄图片数，可以得到以下结论：
- 单张标定板图像： 讨论计算单应性矩阵 H H H需要的特征点数
- 多张标定板图像： 讨论多少个不同的单应性矩阵 H H H（即所需拍摄图像数），可以构建解算内参的方程组
1.所需特征点数
首先来看解算
H
H
H所需的特征点数
H
=
[
h
1
h
2
h
1
]
=
M
⋅
[
r
1
r
2
t
]
=
[
f
x
0
c
x
0
f
y
c
y
0
0
1
]
[
r
1
r
2
t
]
H=\begin{bmatrix} \mathbf{h}_{1} & \mathbf{h}_{2} & \mathbf{h}_{1} \end{bmatrix}=M\cdot\begin{bmatrix} \mathbf{r}_{1} & \mathbf{r}_{2} & \mathbf{t} \end{bmatrix} =\left[\begin{array}{ccc} f_{x} & 0 & c_{x}\\ 0 &f_{y}& c_{y}\\ 0 & 0 & 1 \end{array}\right]\begin{bmatrix} \mathbf{r}_{1} & \mathbf{r}_{2} & \mathbf{t} \end{bmatrix}
H=[h1h2h1]=M⋅[r1r2t]=⎣⎡fx000fy0cxcy1⎦⎤[r1r2t]
H
H
H是一个3x3的矩阵，但由于其中内参矩阵是以齐次坐标的形式进行的表达，因此
H
H
H当中实际只有8个未知量待求解，图像是一个二维的点，一个点包含了
(
u
,
v
)
(u,v)
(u,v)两个方向的坐标值，即可以输出2个方程。因此如果需要求解一副图像的所对应的当前标定板到图像像素的坐标映射关系
H
H
H，至少需要有四个特征点（8个方程）被提取到。当然，实际标定当中还是希望有更多的特征点，可以有更多的方程从而构建超定方程组，冗余帮助标定的的鲁棒性更强，但是注意它们最多也只能提取到8个参数，而一个单应性矩阵
H
H
H实际是有6个外参+4个内参组合而成（10个未知量），因此一副图像即使得到特征点数再多也无法解算得到相机的内参。
在OpenCV当中有封装好的函数来通过非线性最小二乘法计算得到单张标定板图像所对应的单应性矩阵。
Mat cv::findHomography ( InputArray srcPoints,
InputArray dstPoints,
int method = 0,
double ransacReprojThreshold = 3,
OutputArray mask = noArray(),
const int maxIters = 2000,
const double confidence = 0.995
)
2.所需拍摄的标定板图像数
如果只想知道结果，那可以明确的说 ，使用张氏标定至少需要的标定板图像为3张，我这边实际测试的结果8张以上的标定板图像结果和误差才会趋于稳定，因此建议实际标定时使用10张及以上的标定板图像进行标定。
但在这里我还是想从数学推导的角度说明一下，为什么至少需要3张标定图像？ 这部分内容在《学习OpenCV 3》1中也可以找到。首先我们在上一节中提到了，计算一个单应性矩阵 H H H至少需要提供4个特征点（8个方程）来计算得到 H H H中的8个参数，而一个 H H H是由6个外部参数+4个内部参数组合而成。每多一副标定板图像会引入6个新的外部参数，因此很容易得到对于每次参与计算 H H H的8个方程，实际只提供了两个约束用于求解那4个固有的内部参数。那我们来看如何从单应性矩阵 H H H求解得到4个内参
对于 H [ h 1 h 2 h 1 ] = M ⋅ [ r 1 r 2 t ] H\begin{bmatrix} \mathbf{h}_{1} & \mathbf{h}_{2} & \mathbf{h}_{1} \end{bmatrix}=M\cdot\begin{bmatrix} \mathbf{r}_{1} & \mathbf{r}_{2} & \mathbf{t} \end{bmatrix} H[h1h2h1]=M⋅[r1r2t]，根据上一篇笔记的内容，我们可以知道外参旋转矩阵中的列向量相互正交且模相等，即：
（王婆卖瓜）这里用了范数的符号，如果你对范数有兴趣了解的话，可以阅读我的这篇博客【数值分析】学习笔记1——范数与条件数
{
r
1
T
⋅
r
2
=
0
∥
r
1
∥
=
∥
r
1
∥
\left\{\begin{matrix} \mathbf{r}_{1} ^{T}\cdot \mathbf{r}_{2}=0\\\\ \left \| \mathbf{r}_{1} \right \| =\left \| \mathbf{r}_{1} \right \| \end{matrix}\right.
⎩⎨⎧r1T⋅r2=0∥r1∥=∥r1∥根据等式以及
H
H
H的表达式，我们可以得到两个约束：
{
h
1
T
⋅
M
−
T
M
−
1
⋅
h
2
=
0
h
1
T
⋅
M
−
T
M
−
1
⋅
h
1
=
h
2
T
⋅
M
−
T
M
−
1
⋅
h
2
\left\{\begin{matrix} \mathbf{h}_{1} ^{T}\cdot M^{-T}M^{-1}\cdot\mathbf{h}_{2}=0\\\\ \mathbf{h}_{1} ^{T}\cdot M^{-T}M^{-1}\cdot\mathbf{h}_{1}=\mathbf{h}_{2} ^{T}\cdot M^{-T}M^{-1}\cdot\mathbf{h}_{2} \end{matrix}\right.
⎩⎨⎧h1T⋅M−TM−1⋅h2=0h1T⋅M−TM−1⋅h1=h2T⋅M−TM−1⋅h2设
B
=
M
−
T
M
−
1
B=M^{-T}M^{-1}
B=M−TM−1，我们将其展开可以有：
可以看出
B
B
B是一个3x3的对称矩阵，因此只需要求出其中的6个元素即可，设主对角线上方的这6个元素为向量
b
\mathbf{b}
b，则可以将刚刚约束中的式子转换为如下的形式：
h
i
T
⋅
B
⋅
h
j
=
v
i
,
j
T
⋅
b
=
[
h
i
,
1
h
j
,
1
(
h
i
,
1
h
j
,
2
+
h
i
,
2
h
j
,
1
)
h
i
,
2
h
j
,
2
(
h
i
,
3
h
j
,
1
+
h
i
,
1
h
j
,
3
)
(
h
i
,
3
h
j
,
2
+
h
i
,
2
h
j
,
3
)
h
i
,
3
h
j
,
3
]
⋅
b
\begin{aligned} \mathbf{h}_{i} ^{T}\cdot B\cdot\mathbf{h}_{j}&=\mathbf{v}_{i,j}^{T}\cdot\mathbf{b}\\ &=\left[h_{i, 1} h_{j, 1}\left(h_{i, 1} h_{j, 2}+h_{i, 2} h_{j, 1}\right) h_{i, 2} h_{j, 2}\left(h_{i, 3} h_{j, 1}+h_{i, 1} h_{j, 3}\right)\left(h_{i, 3} h_{j, 2}+h_{i, 2} h_{j, 3}\right) h_{i, 3} h_{j, 3}\right]\cdot\mathbf{b} \end{aligned}
hiT⋅B⋅hj=vi,jT⋅b=[hi,1hj,1(hi,1hj,2+hi,2hj,1)hi,2hj,2(hi,3hj,1+hi,1hj,3)(hi,3hj,2+hi,2hj,3)hi,3hj,3]⋅b从而两个约束也可以写成：
[
v
1
,
2
T
(
v
1
,
1
−
v
2
,
2
)
T
]
⋅
b
=
0
\begin{bmatrix} \mathbf{v}_{1,2}^{T}\\ (\mathbf{v}_{1,1}-\mathbf{v}_{2,2})^{T} \end{bmatrix}\cdot \mathbf{b}=0
[v1,2T(v1,1−v2,2)T]⋅b=0即一个包含两个方程的线性方程组，而每一幅图像可以提供一个新的这样的方程组，因此至少需要3张标定板图像（6个方程），才可以求解出
b
\mathbf{b}
b中的6个未知量，从而计算得到4个摄像机系统的内参！
对于求解相机的内参，在OpenCV中同样有有封装好的函数
double cv::calibrateCamera ( InputArrayOfArrays objectPoints,
InputArrayOfArrays imagePoints,
Size imageSize,
InputOutputArray cameraMatrix,
InputOutputArray distCoeffs,
OutputArrayOfArrays rvecs,
OutputArrayOfArrays tvecs,
OutputArray stdDeviationsIntrinsics,
OutputArray stdDeviationsExtrinsics,
OutputArray perViewErrors,
int flags = 0,
TermCriteria criteria = TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 30, DBL_EPSILON)
)
三.标定中的其他注意事项
上面的内容从数学的角度介绍了标定是如何求解得到相机系统针孔成像模型中的各项参数，而对于畸变模型的求解在此就不过多赘述了，可以参照"close range camera calibration"4这篇文章，这也是OpenCV的标定函数所采用的畸变参数求解算法。但这里还是要提几个关键点
- 标定板放置的位置尽可能遍布整个图像（只有特征点覆盖整个视场，畸变校正的鲁棒性才会好）
- 标定时尽量使用辅助光源确保标定平面亮度均匀，过暗和过曝都会导致对于特征点位置提取的不准确
- 标定板的加工精度（特征点间距）和表面平整度对标定有很大的影响，我自制标定板和花钱制作的标定板，重投影误差相差三倍以上
- 尽量保证标定板的姿态多一点，不要总是在一个平面上移动（否则产生的单应性矩阵容易耦合造成标定结果鲁棒性较差）
题外话： 关于手眼标定能否和相机内参标定一起进行的问题
思考： 从算术推导上来看，手眼标定和相机标定是可以同时在一次标定中完成的，但我接触的很多实例中无论眼在手上还是眼在手外，都会偏向于将两者分开来做，我认为最主要的原因还是，如果要保证相机标定的质量好，机械臂带动相机还是标定板，其能运动的范围都是比较有限的，很多关节位置去不到，这样会导致手眼标定时候从相机坐标系到机械臂坐标系的外参鲁棒性较差，因此分开来做也算是一种比较好的选择。
参考引用
凯勒, 布拉德斯基, 刘昌祥, 吴雨培, and 王成龙. 学习OpenCV 3 中文版. 北京: 清华大学出版社, 2018. ↩︎ ↩︎
Zhang Z. A flexible new technique for camera calibration[J]. IEEE Transactions on pattern analysis and machine intelligence, 2000, 22(11): 1330-1334. ↩︎
Duane C B. Close-range camera calibration[J]. Photogramm. Eng, 1971, 37(8): 855-866. ↩︎
标签：begin,end,cdot,标定,相机,mathbf,array,三维重建 来源： https://blog.csdn.net/weixin_43560489/article/details/123170605