- 【一】欧式空间、欧式变换
- 【二】[详细]针孔相机模型、相机镜头畸变模型、相机标定与OpenCV实现
- 【三】仿射变换、投影变换的矩阵形式和特点归纳
- 【四】相机标定
- 【五】边缘检测算子
- 【六】SVD分解
- 【七】GMS算法
- 【八】双边滤波
张氏标定法的原理
1.单应性矩阵H的计算
根据针孔相机模型，可以得到如下表达式：
s
[
u
v
1
]
=
A
[
R
t
]
[
X
W
Y
W
Z
W
1
]
=
A
[
r
1
r
2
r
3
t
]
[
X
W
Y
W
Z
W
1
]
s\left[\begin{array}{l} u \\ v \\ 1 \end{array}\right]=A\left[\begin{array}{ll} R & t \end{array}\right]\left[\begin{array}{c} X_{W} \\ Y_{W} \\ Z_{W} \\ 1 \end{array}\right]=A\left[\begin{array}{llll} r_{1} & r_{2} & r_{3} & t \end{array}\right]\left[\begin{array}{c} X_{W} \\ Y_{W} \\ Z_{W} \\ 1 \end{array}\right]
s???uv1????=A[R?t?]?????XW?YW?ZW?1??????=A[r1??r2??r3??t?]?????XW?YW?ZW?1??????
假设标定板所在的平面为世界坐标系所在的平面，即: Z w = 0 Z_{w}=0 Zw?=0
则上式可以改写为：
s
[
u
v
1
]
=
A
[
R
t
]
[
X
W
Y
W
0
1
]
=
A
[
r
1
r
2
t
]
[
X
W
Y
W
1
]
s\left[\begin{array}{l} u \\ v \\ 1 \end{array}\right]=A\left[\begin{array}{ll} R & t \end{array}\right]\left[\begin{array}{c} X_{W} \\ Y_{W} \\ 0 \\ 1 \end{array}\right]=A\left[\begin{array}{lll} r_{1} & r_{2} & t \end{array}\right]\left[\begin{array}{c} X_{W} \\ Y_{W} \\ 1 \end{array}\right]
s???uv1????=A[R?t?]?????XW?YW?01??????=A[r1??r2??t?]???XW?YW?1????
其中，矩阵H为：
H = A [ r 1 r 2 t ] = [ h 1 h 2 h 3 ] = [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 1 ] H=A\left[\begin{array}{lll} r_{1} & r_{2} & t \end{array}\right]=\left[\begin{array}{lll} h_{1} & h_{2} & h_{3} \end{array}\right]=\left[\begin{array}{llc} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & 1 \end{array}\right] H=A[r1??r2??t?]=[h1??h2??h3??]=???h11?h21?h31??h12?h22?h32??h13?h23?1????
则上式可以展开为：
{
s
u
=
h
11
X
+
h
12
Y
+
h
13
s
v
=
h
21
X
+
h
22
Y
+
h
23
s
=
h
31
X
+
h
32
Y
+
1
\left\{\begin{array}{l} s u=h_{11} X+h_{12} Y+h_{13} \\ s v=h_{21} X+h_{22} Y+h_{23} \\ s=h_{31} X+h_{32} Y+1 \end{array}\right.
????su=h11?X+h12?Y+h13?sv=h21?X+h22?Y+h23?s=h31?X+h32?Y+1?
将(3)中的s带入(1)(2)中得到：
{
(
h
31
X
+
h
32
Y
+
1
)
u
=
h
11
X
+
h
12
Y
+
h
13
(
h
31
X
+
h
32
Y
+
1
)
v
=
h
21
X
+
h
22
Y
+
h
23
\left\{\begin{array}{l} \left(h_{31} X+h_{32} Y+1\right) u=h_{11} X+h_{12} Y+h_{13} \\ \left(h_{31} X+h_{32} Y+1\right) v=h_{21} X+h_{22} Y+h_{23} \end{array}\right.
{(h31?X+h32?Y+1)u=h11?X+h12?Y+h13?(h31?X+h32?Y+1)v=h21?X+h22?Y+h23??
令:
h
′
=
[
h
11
h
12
h
13
h
21
h
22
h
23
h
31
h
32
]
h^{\prime}=\left[\begin{array}{llllllll} h_{11} & h_{12} & h_{13} & h_{21} & h_{22} & h_{23} & h_{31} & h_{32} \end{array}\right]
h′=[h11??h12??h13??h21??h22??h23??h31??h32??]
则改写为矩阵形式为：
[
X
Y
1
0
0
0
?
u
X
?
u
Y
?
u
0
0
0
X
Y
1
?
v
X
?
v
Y
?
v
]
h
′
=
0
\left[\begin{array}{lllllllll} X & Y & 1 & 0 & 0 & 0 & -u X & -u Y & -u \\ 0 & 0 & 0 & X & Y & 1 & -v X & -v Y & -v \end{array}\right] h^{\prime}=0
[X0?Y0?10?0X?0Y?01??uX?vX??uY?vY??u?v?]h′=0
上式可以看作：
S
h
′
=
0
S h^{\prime}=0
Sh′=0
那么矩阵
S
T
S
S^{T}S
STS最小特征值对应的特征向量就是该方程的最小二乘解，再将解归一化得到所需的
h
′
h^{\prime}
h′，从而可以求得
H
H
H。但是由于线性解法所得到的解一般不是最优的解，所以可以选择上面两个等式中的一个，构建评价函数，利用LM算法计算出更高的精度解。
2.相机内外参求解
由于求得的H可能和真实的H相差一个比例因子，因此得到：
[
h
1
h
2
h
3
]
=
λ
A
[
r
1
r
2
t
]
\left[\begin{array}{lll} h_{1} & h_{2} & h_{3} \end{array}\right]=\lambda A\left[\begin{array}{lll} r_{1} & r_{2} & t \end{array}\right]
[h1??h2??h3??]=λA[r1??r2??t?]
其中，内参A矩阵为：
A
=
[
a
γ
u
0
0
β
v
0
0
0
1
]
A=\left[\begin{array}{ccc} a & \gamma & u_{0} \\ 0 & \beta & v_{0} \\ 0 & 0 & 1 \end{array}\right]
A=???a00?γβ0?u0?v0?1????
补充：三阶上三角、下三角矩阵求逆公式：
如果有可逆下三角矩阵：
A
=
[
m
0
0
a
n
0
b
c
h
]
A=\left[\begin{array}{lll} m & 0 & 0 \\ a & n & 0 \\ b & c & h \end{array}\right]
A=???mab?0nc?00h????
则A的逆矩阵为：
A
?
1
=
[
1
m
0
0
?
a
m
n
1
n
0
a
c
?
b
n
m
n
h
?
c
n
h
1
h
]
A^{-1}=\left[\begin{array}{ccc} \frac{1}{m} & 0 & 0 \\ -\frac{a}{m n} &\frac{1}{n} & 0 \\ \frac{a c-b n}{m n h} & -\frac{c}{n h} & \frac{1}{h} \end{array}\right]
A?1=???m1??mna?mnhac?bn??0n1??nhc??00h1?????
如果有可逆上三角矩阵：
B
=
[
m
a
b
0
n
c
0
0
h
]
B=\left[\begin{array}{lll} m & a & b \\ 0 & n & c \\ 0 & 0 & h \end{array}\right]
B=???m00?an0?bch????
则B的逆矩阵为：
B
?
1
=
[
1
m
?
a
m
n
a
c
?
b
n
m
n
h
0
n
ˉ
?
c
n
h
0
0
1
h
]
B^{-1}=\left[\begin{array}{ccc} \frac{1}{m} & -\frac{a}{m n} & \frac{a c-b n}{m n h} \\ 0 & \bar{n} & -\frac{c}{n h} \\ 0 & 0 & \frac{1}{h} \end{array}\right]
B?1=???m1?00??mna?nˉ0?mnhac?bn??nhc?h1?????
那么可以得到内参A矩阵的逆矩阵：
A
?
1
=
[
1
a
?
γ
a
β
γ
v
0
?
u
0
β
a
β
0
1
β
?
v
0
β
0
0
γ
]
A^{-1}=\left[\begin{array}{ccc} \frac{1}{a} & -\frac{\gamma}{a \beta} & \frac{\gamma v_{0}-u_{0} \beta}{a \beta} \\ 0 & \frac{1}{\beta} & -\frac{v_{0}}{\beta} \\ 0 & 0 & \gamma \end{array}\right]
A?1=???a1?00??aβγ?β1?0?aβγv0??u0?β??βv0??γ????
r
1
r_{1}
r1?和
r
2
r_{2}
r2?为单位正交向量，有
r
1
T
r
1
=
r
2
T
r
2
=
1
r_{1}^{T} r_{1}=r_{2}^{T} r_{2}=1
r1T?r1?=r2T?r2?=1，所以上式可以得到两个约束条件：
{
h
1
T
A
?
T
A
?
1
h
2
=
0
h
1
T
A
?
T
A
?
1
h
1
=
h
2
T
A
?
T
A
?
1
h
2
\left\{\begin{array}{l} h_{1}^{T} A^{-T} A^{-1} h_{2}=0 \\ h_{1}^{T} A^{-T} A^{-1} h_{1}=h_{2}^{T} A^{-T} A^{-1} h_{2} \end{array}\right.
{h1T?A?TA?1h2?=0h1T?A?TA?1h1?=h2T?A?TA?1h2??
令：
B
=
A
?
T
A
?
1
=
[
B
11
B
12
B
13
B
21
B
22
B
23
B
31
B
32
B
33
]
=
[
1
a
2
?
γ
a
2
β
v
0
γ
?
u
0
β
a
2
β
?
γ
a
2
β
γ
a
2
β
2
+
1
β
2
?
γ
(
v
0
γ
?
u
0
β
)
a
2
β
2
?
v
0
β
2
v
0
γ
?
u
0
β
a
2
β
?
γ
(
v
0
γ
?
u
0
β
)
a
2
β
2
?
v
0
β
2
?
(
v
0
γ
?
u
0
β
)
2
a
2
β
+
v
0
2
β
2
+
1
]
\begin{array}{l} B=A^{-T} A^{-1}=\left[\begin{array}{lll} B_{11} & B_{12} & B_{13} \\ B_{21} & B_{22} & B_{23} \\ B_{31} & B_{32} & B_{33} \end{array}\right] =\\ \left[\begin{array}{ccc} \frac{1}{a^{2}} & -\frac{\gamma}{a^{2} \beta} & \frac{v_{0} \gamma-u_{0} \beta}{a^{2} \beta} \\ -\frac{\gamma}{a^{2} \beta} & \frac{\gamma}{a^{2} \beta^{2}}+\frac{1}{\beta^{2}} & -\frac{\gamma\left(v_{0} \gamma-u_{0} \beta\right)}{a^{2} \beta^{2}}-\frac{v_{0}}{\beta^{2}} \\ \frac{v_{0} \gamma-u_{0} \beta}{a^{2} \beta} & -\frac{\gamma\left(v_{0} \gamma-u_{0} \beta\right)}{a^{2} \beta^{2}}-\frac{v_{0}}{\beta^{2}} & -\frac{\left(v_{0} \gamma-u_{0} \beta\right)^{2}}{a^{2} \beta}+\frac{v_{0}^{2}}{\beta^{2}}+1 \end{array}\right] \end{array}
B=A?TA?1=???B11?B21?B31??B12?B22?B32??B13?B23?B33?????=????a21??a2βγ?a2βv0?γ?u0?β???a2βγ?a2β2γ?+β21??a2β2γ(v0?γ?u0?β)??β2v0???a2βv0?γ?u0?β??a2β2γ(v0?γ?u0?β)??β2v0???a2β(v0?γ?u0?β)2?+β2v02??+1??????
从上可以看出B是一个对称矩阵，可以用6维向量定义：
b
=
[
B
11
B
12
B
22
B
13
B
23
B
33
]
?
b=\left[\begin{array}{llllll} B_{11} & B_{12} & B_{22} & B_{13} & B_{23} & B_{33} \end{array}\right]^{\top}
b=[B11??B12??B22??B13??B23??B33??]?
假设H的第i列向量可以表示为
h
i
=
[
h
i
1
h
i
2
h
i
3
]
h_{i}=\left[\begin{array}{l} h_{i 1} \\ h_{i 2} \\ h_{i 3} \end{array}\right]
hi?=???hi1?hi2?hi3?????
那么：
h
i
T
B
h
i
=
V
i
j
?
b
j
h_{i}^{T} B h_{i}=V_{i j}^{\top} b_{j}
hiT?Bhi?=Vij??bj?
其中：
V
i
j
=
[
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
]
?
\begin{array}{l} V_{i j} =\left[\begin{array}{lllll} 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} \end{array}\right]^{\top} \end{array}
Vij?=[hi1?hj1??hi1?hj2?+hi2?hj1??hi2?hj2??hi3?hj1?+hi1?hj3??hi3?hj2?+hi2?hj3?+hi3?hj3??]??
则可以得到：
[
V
12
T
V
11
T
?
V
22
T
]
b
=
0
\left[\begin{array}{c} V_{12}^{T} \\ V_{11}^{T}-V_{22}^{T} \end{array}\right] b=0
[V12T?V11T??V22T??]b=0
其中，一个H构成两个约束，因此至少需要三个方程才可以求解出b，从而得到5个内参数：
{
v
0
=
(
B
12
B
13
?
B
11
B
23
)
/
(
B
11
B
22
?
B
12
2
)
λ
=
B
33
?
[
B
13
2
+
v
0
(
B
12
B
13
?
B
11
B
23
)
]
/
B
11
f
u
=
λ
/
B
11
f
v
=
λ
B
11
/
(
B
11
B
22
?
B
12
2
)
s
=
?
B
12
f
u
2
f
v
/
λ
u
0
=
s
v
0
/
f
v
?
B
13
f
u
2
/
λ
\left\{\begin{array}{l} v_{0}=\left(B_{12} B_{13}-B_{11} B_{23}\right) /\left(B_{11} B_{22}-B_{12}^{2}\right) \\ \lambda=B_{33}-\left[B_{13}^{2}+v_{0}\left(B_{12} B_{13}-B_{11} B_{23}\right)\right] / B_{11} \\ f_{u}=\sqrt{\lambda / B_{11}} \\ f_{v}=\sqrt{\lambda B_{11} /\left(B_{11} B_{22}-B_{12}^{2}\right)} \\ s=-B_{12} f_{u}^{2} f_{v} / \lambda \\ u_{0}=s v_{0} / f_{v}-B_{13} f_{u}^{2} / \lambda \end{array}\right.
????????????????v0?=(B12?B13??B11?B23?)/(B11?B22??B122?)λ=B33??[B132?+v0?(B12?B13??B11?B23?)]/B11?fu?=λ/B11??fv?=λB11?/(B11?B22??B122?)?s=?B12?fu2?fv?/λu0?=sv0?/fv??B13?fu2?/λ?
再根据单应性矩阵H和内参矩阵A，利用如下公式，计算每幅图像的外参：
{
r
1
=
λ
A
?
1
h
1
,
r
2
=
λ
A
?
1
h
2
,
r
3
=
r
1
×
r
2
t
=
λ
A
?
1
h
3
,
λ
=
1
∥
A
?
1
h
1
∥
=
1
∥
A
?
1
h
2
∥
\left\{\begin{array}{l} r_{1}=\lambda A^{-1} h_{1}, r_{2}=\lambda A^{-1} h_{2}, r_{3}=r_{1} \times r_{2} \\ t=\lambda A^{-1} h_{3}, \lambda=\frac{1}{\left\|A^{-1} h_{1}\right\|}=\frac{1}{\left\|A^{-1} h_{2}\right\|} \end{array}\right.
{r1?=λA?1h1?,r2?=λA?1h2?,r3?=r1?×r2?t=λA?1h3?,λ=∥A?1h1?∥1?=∥A?1h2?∥1??
由于图像中存在一些噪声，所以矩阵R事实上并不满足正交性质，所以根据最小距离准则取最佳的R解。
标定方法：
| 标定方法 | 优点 | 缺点 | 常见方法 |
|---|---|---|---|
| 标定物标定法 | 可使用任意的相机模型、精度高 | 需要标定物、算法复杂 | Tsai两步法、张氏标定法 |
| 相机自标定法 | 灵活性强、可在线标定 | 精度低、鲁棒性差 | 分层逐步标定、基于Kruppa方程 |
| 主动视觉相机标定法 | 不需要标定物、算法简单、鲁棒性高 | 成本高、设备昂贵、对于运动参数无法预知的情况不适用 | 主动系统控制相机做特定运动 |
| 标定方法 | 相机模型 | 畸变模型 | skewness | 原理 | 优点 | 缺点 |
|---|---|---|---|---|---|---|
| Tsai | 针孔相机模型 | 二阶 | 0 | 给定至少7组特征点的像素坐标和世界坐标，基于径向对齐构建超定线性方程组 | 精度高 | 需要精确的3D测量、耗时、采集数据易受噪声影响、大多数情况不宜使用 |
| zhang | 针孔相机模型 | 二阶+四阶 | 变量 | 计算每个棋盘格的单应性矩阵，至少3个视角的单应性矩阵已知，相机参数可以通过对这些矩阵的线性等式求解得到 | 需要的设备简单（标定板），可以实现较高的精度、灵活、适应性强 | 对噪声敏感，但是可以通过包含更多格子的棋盘格提高精度 |
- 注意：《Requirements for Camera Calibration: Must Accuracy Come with a High Price?》
- 零偏移假设，即s=0，一般情况下是合理的。
- 二阶足以模拟径向畸变模型，四阶在低噪声情况下可取，六阶将降低标定的性能。
- 添加切向畸变，一般可以提高标定的精度。
- 在拍摄时，标定板的图像应尽量覆盖整个视野的1/2左右，并且标定图片一般以20张为宜，实验证明，图片太多会导致参数优化的结果变差
标定板
精度：圆形标定板 > 角点标定板
圆形标定板
- halcon中圆形标定板
- 自制圆形标定板
注意：在中间设置了不对称（既不是轴对称，也不是中心对称）的五个大圆，有两个目的： - 无论怎么摆设标定板，都不会改变圆的排序序号，即在相机标定时，在标定板上建立世界坐标系是唯一的，这在一定程度会提升相机标定的精度。
- 在标定时，由于相机视角原因，即使没有将标定板拍摄完全，也可以根据中间五个大圆，将图像上的圆正确排序。
角点标定板
可以满足VSLAM相机标定的需要
-
棋盘格标定板
-
二维标识码（Apriltag）标定板
