前言
前面的几篇博客中介绍了有关相机标定的基础知识(三维重建学习(1):基础知识:旋转矩阵与旋转向量、三维重建学习(2):相机标定基础)。这次介绍一个十分经典的单目相机标定方法——张正友标定,并给出数学理论推导。
基本方程
模型
我们首先约定如下表示:
二维点坐标:m=[uv],三维点坐标:M=⎡⎣⎢XYZ⎤⎦⎥
将上面的二维点和三维点的坐标表示成齐次形式:
二维点坐标:m~=⎡⎣⎢uv1⎤⎦⎥,,三维点坐标:M~=⎡⎣⎢⎢⎢XYZ1⎤⎦⎥⎥⎥
参考前面的博客:三维重建学习(2):相机标定基础,我们可以建立一个常见的相机小孔成像模型:
s⋅m~=A[Rt]M~
其中:
A=⎡⎣⎢α00cβ0u0v01⎤⎦⎥
s为任意比例因子,
模型平面与图像之间的单应性关系
设
代入原方程中:
s⋅⎡⎣⎢uv1⎤⎦⎥=A[r1r2r3t]⎡⎣⎢⎢⎢XYZ1⎤⎦⎥⎥⎥
假设模型平面在世界坐标系中Z坐标为0,那么
s⋅⎡⎣⎢uv1⎤⎦⎥=A[r1r2r3t]⎡⎣⎢⎢⎢XY01⎤⎦⎥⎥⎥
乘出来都还是0,省略掉这一部分:
s⋅[uv]=A[r1r2t]⎡⎣⎢XY1⎤⎦⎥
此时坐标表示也要变换一下:
二维点坐标:m~=[uv],,三维点坐标:M~=⎡⎣⎢XY1⎤⎦⎥
因此点M和它在图像上的映射点
其中:
H=A[r1r2t]
很显然,H是一个
内参的约束条件
令H=[h1h2h3],h1、h2、h3都是3×1的矩阵,各自对应一列。
则有[h1h2h3]=λA[r1r2t],式中λ是任意的标量。
我们知道旋转矩阵的每一列两两正交,即r1与r2正交。这里不对基础知识赘述,如果有疑问,请查看:旋转矩阵(Rotate Matrix)的性质分析。(吐槽一下:旋转矩阵真的是一个完美的矩阵)
根据r1与r2正交,我们能得到条件:
rT1r2=0,∥r1∥=∥r2∥=1
上面这个条件先放在这里,后面再用。
由前面公式:[h1h2h3]=λA[r1r2t]一一对应得到方程组:
⎧⎩⎨h1=λAr1h2=λAr2h3=λAt
推出:
{r1=λ−1A−1h1r2=λ−1A−1h2
代入前面的条件:
rT1r2=0,∥r1∥=∥r2∥=1
得到:(λ是常数)
rT1r2=λ−ThT1A−TA−1h2λ−1=λ−2hT1A−TA−1h2=0
{∥r1∥2=λ−2hT1A−TA−1h1=1∥r2∥2=λ−2hT2A−TA−1h2=1⇒λ−2hT1A−TA−1h1=λ−2hT2A−TA−1h2
由上面的式子我们可以发现,对于一个给定的单应性矩阵H,对于内参有2个基本的约束条件:
对于H矩阵来说,它是一个
解决相机标定
利用约束条件求出内参矩阵A
好的,还是看一下前面给出的约束条件,注意到中间都有这么一个东西:
我们试着将它表示出来:令B=A−TA−1。
我们在最前面给出了A的定义:
计算A的逆矩阵,步骤就省略了,很简单。得到:
A=⎡⎣⎢⎢⎢⎢⎢1α00−cαβ1β0v0c−u0βαβ−v0β1⎤⎦⎥⎥⎥⎥⎥
计算出B来:
不难发现B是对称的,我们可以使用6个变量来表示出
b=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢B11B12B22B13B23B33⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥
B可以表示为:
与前面表示旋转矩阵时类似,我们假设H的第
对于hTiBhj,我们把前面的公式代入看看:
hTiBhj=[hi1hi2hi3]⋅⎡⎣⎢B11B12B13B12B22B23B13B23B33⎤⎦⎥⋅⎡⎣⎢hj1hj2hj3⎤⎦⎥=(hi1B11+hi2B12+hi3B13)hj1+(hi1B12+hi2B22+hi3B23)hj2+(hi3B13+hi2B23+hi3B33)hj3=hi1hj1B11+(hi2hj1+hi1hj2)B12+hi2hj2B22+(hi3hj1+hi1hj3)B13+(hi3hj2+hi2hj3)B23+hi3hj3B33
我们前面已经给出了一个6维向量b:
那么上面那一堆东西可以表示为:
hTiBhj=vTijb
其中:
vij=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hi1hj1(hi2hj1+hi1hj2)hi2hj2(hi3hj1+hi1hj3)(hi3hj2+hi2hj3)hi3hj3⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥
回到我们前面推导出的约束条件:
{λ−2hT1A−TA−1h2=0λ−2hT1A−TA−1h1=λ−2hT2A−TA−1h2
这两个约束条件可以改写为齐次形式:
{vT12b=0(v11−v22)Tb=0
我们用一个新的矩阵V来表示这两个式子:
这里的V是一个
如果观察了n张图片,那么可以得到n个方程
如果我们求出了唯一解b,那么就可以得到
下面直接给出结果:(已知
利用内参矩阵A求解外参矩阵
通过前面的方法,假设我们已经求到了内参矩阵
我们使用下面的式子表示点M和它在图像上的映射点
s⋅m~=H⋅M~
如果我们已知图片中的点m的坐标,以及点
对于已知的单应性矩阵
[h1h2h3]=λA[r1r2t]
解出外参:
⎧⎩⎨⎪⎪r1=λ−1A−1h1r2=λ−1A−1h2t=λ−1A−1h3
由旋转矩阵的性质得到:
{r3=r1×r2∥r1∥=∥r2∥=1
这些东西整合一下:
⎧⎩⎨⎪⎪⎪⎪⎪⎪r1=λ−1A−1h1r2=λ−1A−1h2r3=r1×r2t=λ−1A−1h3∥λ−1A−1h1∥=∥λ−1A−1h2∥=1⇒λ=∥A−1h1∥=∥A−1h2∥
这样就求出了外参的各项参数。
极大似然估计
普通形式
张正友在论文中提到,前面的这些数学原理和推导并没有太多的物理意义,仅仅是为后面的极大似然优化提供了一个初值。
假设我们得到了模型平面的n幅图片,模型平面上有m个点,假设图像上像素点的噪声服从独立的同一分布,下面给出极大似然优化问题:
∑i=1n∑j=1m∥mij−m^(A,Ri,ti,Mj)∥
其中:m^(A,Ri,ti,Mj)表示的是点Mj在第i幅图像上的投影,旋转矩阵
好了,这是一个非线性优化问题,他衡量的是点的实际坐标与估计值的误差,我们期望它尽可能地小。解决这个优化问题有很多种工具,但不是这里的重点,所以不做赘述。只是有一点需要注意,它需要一个A的初值,我们使用前面那一大堆公式导出一个猜测值赋给它。为什么这么大费周章?很简单,因为这比随机初始化初值收敛的快得多,有助于加快算法。
径向畸变的处理
如前面的博客:三维重建学习(2):相机标定基础中所说,通常畸变有两种:径向畸变、切向畸变。通常切向畸变可以忽略不计,我们只考虑径向畸变。
令
下面给出张正友在论文中给出的公式:
x~=x[1+k1(x2+y2)+k2(x2+y2)2]y~=y[1+k1(x2+y2)+k2(x2+y2)2]
k1、k2是径向畸变参数。
从公式:u~=u0+αx~+cy~、v~=v0+βy~(注:由内参矩阵的参数导出),得到:
u~=u+(u−u0)[k1(x2+y2)+k2(x2+y2)2]v~=v+(v−v0)[k1(x2+y2)+k2(x2+y2)2]
从上面的式子我们可以得到两个方程:
u~−u=(u−u0)[k1(x2+y2)+k2(x2+y2)2]v~−v=(v−v0)[k1(x2+y2)+k2(x2+y2)2]
表示成矩阵形式:
[(u−u0)(x2+y2)(v−v0)(x2+y2)(u−u0)(x2+y2)2(v−v0)(x2+y2)2]⋅[k1k2]=[u~−uv~−v]
n幅图像中各有
用矩阵表示来表示这个方程组:D⋅k=d,其中k=[k1k2]。
使用最小二乘法求解出k,直接套公式:
如此求出了畸变参数后,我们可以回到前面的优化问题中再求解其他参数。
完整形式
下面给出完整的极大似然估计优化问题:
∑i=1n∑j=1m∥mij−m^(A,k1,k2,Ri,ti,Mj)∥
其中:m^(A,ki,kj,Ri,ti,Mj)表示的是点Mj在第i幅图像上的投影。
与前面类似,这里依然回一个非线性优化问题。其中内参矩阵
到这里,整个算法的数学原理已经推导完毕了。
在论文中,张正友还给出了相机标定的程序流程:
个人感觉用中文翻译意思上总是会有点偏差,姑且还是给出翻译后的流程:
建议采用如下校准步骤:
1. 打印一个图案,并把它贴到一个平面上;
2. 通过移动平面或者相机,从不同的方向拍摄一些图片;
3. 检测图片中的特征点;
4. 采用3.1节所述的方法,估计5个内参和全部的外参;
5. 使用最小二乘法(13)估计径向畸变系数;
6. 最小化优化问题(14),改进所有的参数
这里用到的数学公式全部都在前面介绍了,理应很清楚了,不做赘述。下一篇博客再进行程序实现。
参考资料:
[图像]摄像机标定(2) 张正友标定推导详解张正友标定法的真实理解Flexible Camera Calibration By Viewing a Plane From Unknown Orientations(Zhengyou Zhang)