张正友标定数学原理推导
整体概述:
1.先求解单应性矩阵,单应性矩阵是从一个平面 到另一个平面的投影映射
2.通过单应性矩阵求解相机内参参数
3.求解外参参数
4.几何解释
求解单应性矩阵
其实质对应空间坐标系旋转平移变换后通过相机内参矩阵映射到像素坐标。
A称为相机的内参矩阵,内参矩阵取决于相机的内部参数。其中, f f f 是通过欧拉积分为像距, d X , d Y dX,dY dX,dY分别表示 X , Y X,Y X,Y方向上的一个像素在相机感光板上的物理长度(即一个像素在感光板上是多少毫米), u 0 , v 0 u_0,v_0 u0,v0分别表示相机感光板中心在像素坐标系下的坐标, θ \theta θ表示感光板的横边和纵边之间的角度( 90 ° 90° 90°表示无误差)。
我们的关注点不是定义在所有空间的坐标,只是定义在我们所观察的平面的坐标,我们可以选择物体平面为Z=0,
写为单应性矩阵形式
在opencv中提到单应性矩阵中只有8个自由参数,我们进行标准化使H33 = 1(通常是可行的,除了非常罕见的奇异情况H33=0)。单应性矩阵的缩放可以应用于第九个单应性矩阵的参数,但通常倾向于将整个单应性矩阵乘以比例因子来缩。
我个人的理解就是把H33归一化为1,因为前面的比例因子不会对映射关系产生影响,可以将其提到比例因子中。
将单应性矩阵拆分方程求解,将第三行乘上ui与第一行消去比例因子s,第二行与第一行操作相同得到
对应的一组点有两个方程,所以棋盘格四个角点8个方程即可解的8个未知数。
求解相机内参参数
通过两个约束求解,r1和r2是单位正交向量,用h表示
令B为
令b为
B为对称矩阵,所以b的6维向量可以囊括所有B元素
上面两个约束可以表示为通用形式
上面的一幅标定板图像取得的约束等式,假如有n幅图像,则
V
b
=
0
Vb = 0
Vb=0
n幅图即2n×6的矩阵,b有6个未知数用三幅图以上会得到6个以上方程即可求解。
其中,V是一个2n×6的矩阵,bb是一个6维向量,所以
- 当n≥3,可以得到bb的唯一解;
- 当n=2,则可以假设扭曲参数γ=0γ=0作为额外的约束条件
- 当n=1,则值能计算两个相机的内参数
对于方程Vb=0可以使用SVD求得其最小二乘解。对VTV进行SVD分解,其最小特征值对应的特征向量就是Vb=0的最小二乘解,从而求得矩阵B。由于这里得到的B的估计值是在相差一个常量因子下得到的,所以有:
B
=
λ
(
A
−
)
T
A
B=λ(A^-)^TA
B=λ(A−)TA
通过b也就知道了B,通过B可以提取内参参数
上述公式就是内参的解析解,内参参数获得后,外参矩阵便可以求解。
求解外参矩阵
利用上面H与A的关系,R矩阵的三向量为正交向量可得
且r1和r2是单位向量,所以
最大似然估计
上面使用最小二乘法得到估计得到的解,并没有物理上的实际意义,。为了进一步增加标定结果的可靠性,可以使用最大似然估计(Maximum likelihood estimation)来优化上面估计得到的结果。
假设同一相机从n个不同的角度的得到了n幅标定板的图像,每幅图像上有m个像点。Mij表示第ii幅图像上第j个像点对应的标定板上的三维点,则
一般图像噪声默认为高斯噪声,所以使用最大似然估计进行优化,角点mij的概率密度函数为
f
(
m
i
j
)
=
1
2
π
e
−
(
(
^
m
)
(
K
,
R
i
,
t
i
,
M
i
j
)
)
σ
2
f(m_{ij})=\frac{1}{\sqrt{2\pi}}e^{\frac{-(\hat(m)(K, R_i, t_i, M_{ij}))}{\sigma^2}}
f(mij)=2π1eσ2−((^m)(K,Ri,ti,Mij))
构造似然函数
L
(
A
,
R
i
,
t
i
,
M
i
j
)
=
∏
i
=
1
,
j
=
1
n
,
m
f
(
m
i
j
)
=
1
2
π
e
−
∑
i
=
1
n
∑
j
=
1
m
(
m
^
(
K
,
R
i
,
t
i
,
M
i
j
)
−
m
i
j
)
2
σ
2
L(A,R_i,t_i,M_{ij}) = \prod^{n,m}_{i=1,j=1}f(m_{ij})=\frac{1}{\sqrt{2\pi}}e^{\frac{-\sum^n_{i=1}\sum^m_{j=1}(\hat{m}(K,R_i,t_i,M_{ij})-m_{ij})^2}{\sigma^2}}
L(A,Ri,ti,Mij)=i=1,j=1∏n,mf(mij)=2π1eσ2−∑i=1n∑j=1m(m^(K,Ri,ti,Mij)−mij)2
为了能够让L取得最大值,需要最小化下面的值
∑
i
=
1
n
∑
j
=
1
m
∥
m
(
K
,
R
i
,
t
i
,
M
i
j
)
−
m
i
j
∥
2
∑i=1n∑j=1m∥m^(K,Ri,ti,Mij)−mij∥2
∑i=1n∑j=1m∥m(K,Ri,ti,Mij)−mij∥2
这是一个非线性优化问题,可以使用Levenberg-Marquardt的方法,利用上面得到的解作为初始值,迭代得到最优解。
几何解释
在张正友标定法原文中作者给出了几何解释,利用绝对圆锥曲线将两个旋转向量的约束关系联系起来。
绝对圆锥曲线内容参考
两个约束关系
r
1
T
⋅
r
2
=
0
r_1^T\centerdot r_2 = 0
r1T⋅r2=0
r
1
T
⋅
r
1
=
r
2
T
⋅
r
2
r_1^T\centerdot r_1 = r_2^T\centerdot r_2
r1T⋅r1=r2T⋅r2
绝对圆锥曲线参考博客Absolute Conic
三维坐标点表示为齐次坐标为四维向量,平面方程表达式写为上式并用来表示相机坐标系
在齐次坐标系中无穷远点收为同一个点,所以该平面与无穷远平面交于一条线,
[
r
1
0
]
与
[
r
2
0
]
\left[ \begin{matrix} r_1 \\ 0 \end{matrix} \right]与 \left[ \begin{matrix} r_2 \\ 0 \end{matrix} \right]
[r10]与[r20]
是该相交直线上两个特殊的点,所以该直线上任何点可以用这两个点进行表示
x
∞
=
a
[
r
1
0
]
+
b
[
r
2
0
]
=
[
a
r
1
+
b
r
2
0
]
x_\infty = a\left[ \begin{matrix}r_1\\0 \end{matrix}\right]+b\left[ \begin{matrix}r_2\\0 \end{matrix}\right] = \left[\begin{matrix}ar_1+br_2\\0\end{matrix}\right]
x∞=a[r10]+b[r20]=[ar1+br20]
求该直线与绝对圆锥曲线的交点,
x
∞
T
x
∞
=
0
x_\infty^Tx_\infty = 0
x∞Tx∞=0 或
a
2
+
b
2
=
0
a^2+b^2 = 0
a2+b2=0,求解可得
x
∞
=
a
[
r
1
±
i
r
2
0
]
x_\infty = a\left[\begin{matrix}r_1\pm ir_2\\0\end{matrix}\right]
x∞=a[r1±ir20]
这对复共轭点的意义在于它们对欧几里德变换是不变的。 它们在图像平面上的投影(由比例因子决定)为
该点在绝对圆锥曲线图上,描述为
A
−
T
A
−
1
A^{-T}A^{-1}
A−TA−1,得到
这就要求上面的两个约束条件满足,这就是其几何意义。
齐次坐标本来就是为了解决平移矩阵的加法问题,使其能够像旋转矩阵一样对坐标系进行乘法,本来的几何意义就比较抽象,对上文的解释也是自己的部分理解并不透彻,只能从齐次坐标中两平行线无穷远交汇于一点来扩充理解。