|
|
5 基于2D标定物的标定方法
基于2D标定物的标定方法，原理与基于3D标定物相同，只是通过相机对一个平面进行成像，就可得到相机的标定参数，由于标定物为平面，本身所具有的约束条机，相对后者标定更为简单。经典算法为Z. Zhang(PAMI, 2000) A Flexible New Technique for Camera Calibration。其算法已经被收入Opencv(2004)，最常用的标定图案是棋盘格图案，如下图：
5.1 单应性矩阵
对于2D标定平面，抑或称为标定板，不妨假设，平面上点的增广齐次向量[X,Y,Z,1]
[X,Y,Z,1]
满足Z=0
Z=0
，同时对于旋转矩阵R
R
的每一列表示为一个列向量记为ri
ri
，则根据相机的投影方程：
s⎡⎣⎢uv1⎤⎦⎥=A[r1 r2 r3 t]⎡⎣⎢⎢⎢XY01⎤⎦⎥⎥⎥=A[r1 r2 t]⎡⎣⎢XY1⎤⎦⎥(1)s⎡⎣⎢uv1⎤⎦⎥=A[r1 r2 r3 t]⎡⎣⎢⎢⎢XY01⎤⎦⎥⎥⎥=A[r1 r2 t]⎡⎣⎢XY1⎤⎦⎥(1)
因为点[X,Y]T
[X,Y]T
仍表示的是三维空间点坐标，只是由于标定使用平面的特殊性，将Z=0
Z=0
省略，因此我们仍然使用统一的M
M
符号表示，与其对应M~=[X,Y,1]T
M~=[X,Y,1]T
表示该点的其次向量，则公式(1)
(1)
可简记为：
sm~=HM~ with H=A[r1 r2 t](2)sm~=HM~ with H=A[r1 r2 t](2)
H
H
被称为单应性矩阵（Homography matrix）。
5.2 内参约束条件
对于单应性矩阵H
H
，我们也将其按照列向量的方式表示：H=[h1h2h3]
H=[h1 h2 h3]
， 则有：
[h1 h2 h3]=λA[r1 r2 t](3)[h1 h2 h3]=λA[r1 r2 t](3)
其中，λ
λ
是一个任意的系数，因为r1,r2
r1,r2
是正交向量，因此有：
rT1r2=hT1A−TA−1h2=0(4)rT1r2=hT1A−TA−1h2=0(4)
rT1r1=hT1A−TA−1h1=rT2r2=hT2A−TA−1h2(5)rT1r1=hT1A−TA−1h1=rT2r2=hT2A−TA−1h2(5)
对于一个单应性矩阵，公式(4,5)
(4,5)
是两个内参的基本约束。单应性矩阵有9个元素，但可以由8个独立不相关的元素表示，也就是说投影变换有8个自由度，而我们只有6个外参元素（3个旋转元素和3个平移元素），因此两个内参约束条件的意义就在于此。在原理简介（三）中，我们已经了解到A−TA−1
A−TA−1
描述的就是IAC（Image of the absolute conic），后文将对此进行解释。
5.3 几何解释
现在让我们来分析公式(4,5)
(4,5)
与绝对圆锥曲线的关系。在像空间坐标系中，存在下面的等式：
[r3rT3t]T⎡⎣⎢⎢⎢xyzw⎤⎦⎥⎥⎥=0(6)[r3rT3t]T⎡⎣⎢⎢⎢xyzw⎤⎦⎥⎥⎥=0(6)
当w=0
w=0
时就是我们前面解释的改点位于无穷远处。我们想象标定板所在平面与该无穷远处平面相交于一点，则点[r10]
[r10]
和[r20]
[r20]
是相交直线上的两个特殊点，线上的其他点都可以用这两个点线性表示：
x∞=a[r10]+b[r20]=[ar1+br20](7)x∞=a[r10]+b[r20]=[ar1+br20](7)
现在让我们看一下上述的交线和绝对圆锥曲线的交点，原理简介（三）已经介绍点x∞
x∞
满足：xT∞x∞=0
xT∞x∞=0
，也即(ar1+br2)T(ar1+br2)=0⇒a2+b2=0
(ar1+br2)T(ar1+br2)=0 ⇒ a2+b2=0
（r1
r1
与r2
r2
正交）。因此b=±ai,其中i2=−1
b=±ai, 其中i2=−1
，公式(7)
(7)
可写为：
x∞=a[r1±ir20](8)x∞=a[r1±ir20](8)
其在图像中的投影点为：
m∞=A(r1±ir2)=h1±ih2(9)m∞=A(r1±ir2)=h1±ih2(9)
因为点m∞
m∞
位于IAC上，有：
(h1±ih2)TA−TA−1(h1±ih2)=0(10)(h1±ih2)TA−TA−1(h1±ih2)=0(10)
因此等式(10)
(10)
左边的实部与虚部都为0，也就是公式(4,5)
(4,5)
这两条约束条件。
5.4 闭合解
现在开始讲解如何高效求解本方法的相机标定问题。做法是，首先获得分析解，然后初始估计值基于最大似然准则进行非线性优化，这些都将逐步进行讲解。
令：
B=A−TA−1=⎡⎣⎢B11B12B13B12B22B23B13B23B33⎤⎦⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢1α2−γα2βv0γ−u0βα2β−γα2βγ2α2β2+1β2−γ(v0γ−u0β)α2β2−v0β2v0γ−u0βα2β−γ(v0γ−u0β)α2β2−v0β2(v0γ−u0β)2α2β2+v20β2+1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥(11)B=A−TA−1=⎡⎣⎢B11B12B13B12B22B23B13B23B33⎤⎦⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢1α2−γα2βv0γ−u0βα2β−γα2βγ2α2β2+1β2−γ(v0γ−u0β)α2β2−v0β2v0γ−u0βα2β−γ(v0γ−u0β)α2β2−v0β2(v0γ−u0β)2α2β2+v20β2+1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥(11)
因此可以看出，B
B
是一个对称矩阵，由六个元素组成：
b=[B11,B12,B13,B22,B23,B33]T(12)b=[B11,B12,B13,B22,B23,B33]T(12)
对于H=[h1,h2,h3]
H=[h1,h2,h3]
的第i
i
列记为：hi=[hi1,hi2,hi3]T
hi=[hi1,hi2,hi3]T
，则有：
hTiBhj=vTijb(13)hTiBhj=vTijb(13)
其中，vij=[hi1hj1,hi1hj2+hi2hj1,hi2hj2,hi3hj1+hi1hj3,hi3hj2+hi2hj3,hi3hj3]T
vij=[hi1hj1,hi1hj2+hi2hj1,hi2hj2,hi3hj1+hi1hj3,hi3hj2+hi2hj3,hi3hj3]T
。因此，对于公式(4,5)
(4,5)
可以表达为齐次方程：
[vT12(v11−v22)T]b=0(14)[vT12(v11−v22)T]b=0(14)
如果采集了n
n
张图像，可以列出方程组：
Vb=0(15)Vb=0(15)
其中，V
V
是一个2n×6
2n×6
的矩阵。如果n⩾3
n⩾3
，一般就能获得b
b
的唯一解；如果n=2
n=2
，我们可以利用偏斜度γ=0
γ=0
的约束条件，将[0,1,0,0,0,0]b=0
[0,1,0,0,0,0]b=0
添加到方程组(15)
(15)
中；如果n=1
n=1
那么就只能获得相机的内参。
一旦b
b
确认后，就可以获得相机内参。将公式(11)
(11)
中的矩阵B
B
添加一个任意系数：B=λA−TA
B=λA−TA
，可得到：
v0=(B12B13−B11B23)/(B11B22−B212)(16)v0=(B12B13−B11B23)/(B11B22−B212)(16)
λ=B33−[B213+v0(B12B13−B11B23)]/B11(17)λ=B33−[B213+v0(B12B13−B11B23)]/B11(17)
α=λ/B11−−−−−√(18)α=λ/B11−−−−−√(18)
β=λB11/(B11B12−B212)−−−−−−−−−−−−−−−−−√(19)β=λB11/(B11B12−B212)−−−−−−−−−−−−−−−−−√(19)
γ=−B12α2β/λ(20)γ=−B12α2β/λ(20)
u0=γv0/α−B13α2/λ(21)u0=γv0/α−B13α2/λ(21)
内参矩阵获得后，根据公式(2)
(2)
外参也就很快就可以得到：
r1=λ′A−1h1, r2=λ′A−1h2, r3=r1×r2, t=λ′A−1h3(22)r1=λ′A−1h1, r2=λ′A−1h2, r3=r1×r2, t=λ′A−1h3(22)
其中，λ′=1/∥A−1h1∥=1/∥A−1h2∥
λ′=1/∥A−1h1∥=1/∥A−1h2∥
。
由于噪声的存在，这样求解的矩阵R
R
一般并不具备旋转矩阵的特性，更好的求解方法是通过奇异值分解的方法，可以参考Z. Zhang的论文。
5.5 最大似然优化
与基于3D标定物的优化方法类似，这里仍然使用最大似然法进行优化，也就是基于噪声是不相关且独立分布的的假设，对于n
n
张图像，其中包含m
m
个点，优化方程为：
∑i=0n∑j=0m∥mij−m^(A,Ri,ti,Mj)∥2(23)∑i=0n∑j=0m∥mij−m^(A,Ri,ti,Mj)∥2(23)
其中，m^(A,Ri,ti,Mj)
m^(A,Ri,ti,Mj)
是三维点Mj
Mj
在图像i
i
中的投影点。非线性优化过程是，首先需要获得内参矩阵A
A
的估算值，然后利用上面的描述的求解方法，获得{Ri,ti|i=1,…,n}
{Ri,ti|i=1,…,n}
，利用LM迭代优化方法，使得公式(23)
(23)
的求和最小化。
5.6 镜头畸变
在原理简介（四）4.6中已经对镜头畸变模型进行了讲解，这里与之原理一样。以径向畸变为例，为了获得较为理想准确的像点坐标值，则有方程：
[(u−u0)(x2+y2)(v−v0)(x2+y2)(u−u0)(x2+y2)2(v−v0)(x2+y2)2)][k1k2]=[u^−uv^−v](24)[(u−u0)(x2+y2)(v−v0)(x2+y2)(u−u0)(x2+y2)2(v−v0)(x2+y2)2)][k1k2]=[u^−uv^−v](24)
同样当m
m
个点在n
n
张图像中时，可以组成方程组Dk=d
Dk=d
，其中k=[k1k2]
k=[k1k2]
，其线性最小二乘解为：
k=(DTD)−1DTd(25)k=(DTD)−1DTd(25)
当获得(k1,k2)
(k1,k2)
后，我们就可以将公式(23)
(23)
的非线性优化调整为：
∑i=0n∑j=0m∥mij−m^(A,k1,k2,Ri,ti,Mj)∥(26)∑i=0n∑j=0m∥mij−m^(A,k1,k2,Ri,ti,Mj)∥(26)
1 打印出一张标定图并贴到一个平面上；
2 通过移动相机或者标定平面采集不同位置、不同方向的标定板图像；
3 特征点检测；
4 估算内参（公式(16−21)
(16−21)
），然后通过闭合解得到外参（公式(22)(22)）；
5 通过最小二乘线性估算畸变系数；
6 使公式(26)(26)的代数和最小，优化所有参数。 |
|
发表于 2023-1-14 12:33:55