传统的、广角的和鱼眼的镜头的通用模型和标定方法
作者: Juho Kannala and Sami S. Brandt
原文: A Generic Camera Model and Calibration Method for Conventional,
Wide-Angle, and Fish-Eye Lenses
摘要
鱼眼镜头在那种需要非常宽广视角的应用中非常便利,但又由于缺乏准确性,通用性以及易于使用的标定程序,使得它在以测量为目的的应用中受到限制。因此, 我们提出了一个通用的相机模型, 它适用于鱼眼镜头相机以及传统和广角镜头相机, 并给出了一种估计模型参数的校准方法。该标定方法的精确度跟以前最先进的标定方法精度相当。
索引:相机模型,相机标定,镜头畸变,鱼眼镜头,广角镜头
I. 介绍
对于窄角镜头和广角镜头1 2 3来说,带有畸变的针孔模型能够很好的近似,但是对于鱼眼镜头仍然不适合。鱼眼镜头被设计成覆盖相机前方的全部半球形区域,视角非常大,接近1800,更进一步说,通过透视投影投影整个半球形视角到有限的图像平面上是不可行的。因此鱼眼镜头被设计成其他的相机投影模型。这就是为什么鱼眼镜头的固有(内在)畸变不能仅仅视为针孔模型偏差的原因4。
5 6 7已经进行了一些努力来模拟具有不同模型的鱼眼镜片的径向对称畸变。这些方法的思想都是变换鱼眼图像使其遵循针孔模型。畸变模型的参数通过将不是直线的线经过强制变化成直线来估计8 9,该方法没有给出完整的标定,它们仅仅可以 “纠正” 图像以遵循针孔模型, 但当需要知道与图像点相对应的反投影射线的方向时, 其适用性是有限的。10 11意在标定鱼眼镜头。然而, 由于需要激光束或圆柱形校准对象, 这些方法在实践中略显繁琐。
最近, 首次出现了鱼眼镜头相机的自动校准方法12 13 14。Miˇcuˇ s´ık和 B. and Pajdla12提出了极几何(epipolar geometry)和全向相机模型(omnidirectional camera mode)同时线性估计的方法,.Claus和Fitzgibbon14提出了类似的相机运动和镜头几何的同时在线估计方法,Thirthala和Pollefeys13使用径向一维相机的多视图几何来估计非参数相机模型。另外, Barreto 和 Daniilidis15介绍了一种用于校正广角透镜变形的径向基本矩阵。 然而,这些方法的重点是在自动校准技术,而不是在实际镜头的精确建模。
本文主要集中在实际的精准相机模型 2,并提出了一个新颖的标定鱼眼相机的方法,该方法需要相机观测一个平面模板,例如棋盘格。该标定方法基于一个通用摄像机模型,该相机模型适用于不同类型的全向摄像机以及传统摄像机。首先,在第二节中,我们提出了相机模型,在第三节中,通过比较不同的投影模型,从理论上证明了这一点。在第四节中,我们描述了一个估计相机参数的程序,并在第五节和第六节中陈述和讨论了实验结果。
II. 通用相机模型
由于透视投影模型不适合鱼眼镜头,我们使用更加灵活的径向对称投影模型。这个模型在第二节A部分介绍,接着第二节B部分扩展了径向对称到非对称畸变项,回投影(back-projection)的计算在第二节C部分描述。
A. 径向对称模型
针孔相机的透视投影可以使用下面的公式描述:
r
=
f
tan
θ
(
i
.
透
视
投
影
)
(1)
r = f\tan\theta \qquad (i. 透视投影) \tag{1}
r=ftanθ(i.透视投影)(1)
其中
θ
\theta
θ是主轴(principal axis)和入射线的夹角,
r
r
r是图像点和主点之间的距离,
f
f
f是焦距,然而鱼眼镜头通常被设计成服从下面的投影:
r
=
2
f
tan
(
θ
/
2
)
(
i
i
.
立
体
投
影
)
(2)
r = 2f\tan(\theta/2) \qquad (ii. 立体投影) \tag{2}
r=2ftan(θ/2)(ii.立体投影)(2)
r = f θ ( i i i . 等 距 投 影 ) (3) r = f\theta \qquad (iii. 等距投影) \tag{3} r=fθ(iii.等距投影)(3)
r = 2 f sin ( θ / 2 ) ( i v . 等 角 投 影 ) (4) r=2f\sin(\theta/2) \qquad (iv. 等角投影 \tag{4}) r=2fsin(θ/2)(iv.等角投影)(4)
r = f sin ( θ ) ( v . 正 交 投 影 ) (5) r=f\sin(\theta) \qquad (v. 正交投影 \tag{5}) r=fsin(θ)(v.正交投影)(5)
可能最常见的模型是等距投影,不同投影的描述如fig.1(a)所示,针孔相机和鱼眼相机的不同如fig.1(b)所示。
然而实际的镜头不可能是上面提到的精确的设计模型。从自动标定的视角来说,如何我们能够仅仅使用一种模型,它能够适用各种镜头那将非常的有用。因此,我们考虑一般的投影形式:
r ( θ ) = k 1 θ + k 2 θ 3 + k 3 θ 5 + k 4 θ 7 + k 5 θ 9 + . . . (6) r(\theta)=k_1\theta + k_2 \theta^3 + k_3 \theta^5 + k_4 \theta^7 + k_5 \theta^9 + ... \tag{6} r(θ)=k1θ+k2θ3+k3θ5+k4θ7+k5θ9+...(6)
为了不失一般性,去掉了偶次方项。这是因为可以把 r r r作为奇数函数扩展到负向的一侧(因为奇数函数可以张成(span)连续奇数集)。为了计算等式(6),我们需要固定项的数量。我们发现,从 θ \theta θ 5项到9项对于不同的投影曲线都能够很好自由度的近似。因此,我们相机模型径向对称部分包含5个参数, k 1 , k 2 , k 3 , k 4 , k 5 k_1, k_2, k_3, k_4, k_5 k1,k2,k3,k4,k5。
让
F
F
F表示从入射光线到归一化图像坐标的映射:
[
x
y
]
=
r
(
θ
)
[
cos
ψ
sin
ψ
]
=
F
(
Φ
)
(7)
\begin{bmatrix} x & \\ y \end{bmatrix} =r(\theta) \begin{bmatrix} \cos\psi \\ \sin\psi \end{bmatrix} =\mathcal{F(\Phi)} \tag{7}
[xy]=r(θ)[cosψsinψ]=F(Φ)(7)
其中
r
(
θ
)
r(\theta)
r(θ)包含了等式6前5项,
Φ
=
(
θ
,
ψ
)
T
\Phi=(\theta,\psi)^T
Φ=(θ,ψ)T是入射线的方向。对于真实的镜头,参数
k
i
k_i
ki的值使得
r
(
θ
)
r(\theta)
r(θ)在区间
[
0
,
θ
m
a
x
]
[0,\theta_{max}]
[0,θmax]上单调增加,其中
θ
m
a
x
\theta_{max}
θmax表示最大视角。因此,当计算
F
F
F的逆时,我们通过数值找到九阶多项式的根,然后选择0到
θ
m
a
x
之
间
的
实
根
\theta_{max}之间的实根
θmax之间的实根求解
θ
\theta
θ。
B. 完整模型(Full Model)
真实的镜头可能与精确的径向对称偏离,因此,我们提供了一个非对称部分。例如,镜头元件可能不是精确的对准,导致投影不是精确的径向对称。对于传统镜头,这种失真(畸变)叫做偏心(偏离中心)失真。然而,光学系统中还存在其他的缺陷源,其中一些可能难以建模。例如,图像平面可能相对主轴倾斜,或者各个元器件可能不是精确地径向对称。因此,我们没有试图单独模拟光学系统中的各种不同的物理现象,而是提出了灵活的数学失真模型,该模型适合观察结果。
为了获得广泛的应用,灵活的模型,我们提出了两种畸变项,一种畸变项在径向方向上起作用:
δ
r
(
θ
,
ψ
)
=
(
l
1
θ
+
l
2
θ
3
+
l
3
θ
5
)
(
i
1
cos
ψ
+
i
2
sin
ψ
+
i
3
cos
2
ψ
+
i
4
sin
2
ψ
)
(8)
\delta_r(\theta, \psi) = (l_1\theta+l_2\theta^3 +l_3\theta^5)(i_1\cos\psi+i_2\sin\psi+i_3\cos2\psi+i_4\sin2\psi) \tag{8}
δr(θ,ψ)=(l1θ+l2θ3+l3θ5)(i1cosψ+i2sinψ+i3cos2ψ+i4sin2ψ)(8)
一种在切向方向上起作用:
δ
t
(
θ
,
ψ
)
=
(
m
1
θ
+
m
2
θ
3
+
m
3
θ
5
)
(
j
1
cos
ψ
+
j
2
sin
ψ
+
j
3
cos
2
ψ
+
j
4
sin
2
ψ
)
(9)
\delta_t(\theta, \psi)=(m_1\theta + m_2\theta^3+m_3\theta^5)(j_1\cos\psi+j_2\sin\psi+j_3\cos2\psi+j_4\sin2\psi) \tag{9}
δt(θ,ψ)=(m1θ+m2θ3+m3θ5)(j1cosψ+j2sinψ+j3cos2ψ+j4sin2ψ)(9)
畸变函数中的 θ \theta θ变量和 ψ \psi ψ变量是分开的。因为任何连续 2 π 2\pi 2π周期的傅里叶级数收敛于 L 2 L_2 L2范数,并且连续奇函数都可以用一系列奇数多项式表示,原则上可以通过简单的增加更多项来模拟任何类型的连续失真,等式8 和等式9就是那样,都有7个参数模拟。增加径向和切向畸变到等式(7),我们可或得失真(畸变)的坐标 x d = ( x d , y d ) T {\mathbf{x_d}} = (x_d, y_d)^T xd=(xd,yd)T,即:
x
d
=
r
(
θ
)
u
r
(
ψ
)
+
Δ
r
(
θ
,
ψ
)
u
r
(
ψ
)
+
Δ
t
(
θ
,
ψ
)
u
ψ
(
ψ
)
(10)
\mathbf{{x_d}} = r(\theta)\mathbf{u_r}(\psi)+\Delta_r(\theta, \psi)\mathbf{u_r}(\psi)+\Delta_t(\theta, \psi)\mathbf{u_{\psi}}(\psi) \tag{10}
xd=r(θ)ur(ψ)+Δr(θ,ψ)ur(ψ)+Δt(θ,ψ)uψ(ψ)(10)
其中,
u
r
(
ψ
)
\bf{u_r(\psi)}
ur(ψ)和
u
ψ
(
ψ
)
\bf{u_{\psi}}(\psi)
uψ(ψ)是径向和切向方向上的单位向量。为了获取完整的相机模型,我们仍然需要变换传感器平面坐标到图像像素坐标。假设像素坐标系是正交的,我们将从下面的表达式中得到像素坐标
(
u
,
v
)
T
(u, v)^T
(u,v)T
[ u v ] = [ m u 0 0 m v ] [ x d y d ] [ u 0 v 0 ] = A ( x d ) (11) \begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} m_u & 0 \\ 0 & m_v \end{bmatrix} \begin{bmatrix} x_d \\ y_d \end{bmatrix} \begin{bmatrix} u_0 \\ v_0 \end{bmatrix}= \mathcal{A({\mathbf{x_d})}} \tag{11} [uv]=[mu00mv][xdyd][u0v0]=A(xd)(11)
其中 ( u 0 , v 0 ) T (u_0, v_0)^T (u0,v0)T是主点坐标, m u m_u mu和 m v m_v mv分别是水平方向和垂直方向上单位距离的像素的个数。结合等式(10)和等式(11),我们有前向相机模型:
m = P c ( Ψ ) (12) \bf{m} = \mathcal{P_c(\Psi)} \tag{12} m=Pc(Ψ)(12)
其中, m = ( u , v ) T m = (u,v)^T m=(u,v)T,full 相机模型包含23个参数,使用 p 23 p_{23} p23表示。由于模型的不对称部分非常灵活,因此,有时使用更少的的相机模型参数是合理的,避免了过度的拟合。例如,如果控制点(control point)没有覆盖整个图像区域就是这种情况。这样我们去掉不对称部分,给出了 p 9 p_9 p9的相机模型:径向部分对应到等式7,仿射部分对应到等式11的四个参数。
C. Backward Model (反投影模型)
上面我们已经描述了前向相机模型 P c P_c Pc,实际上,我们也需要知道反投影模型:
Ψ
=
P
c
−
1
(
m
)
(13)
\Psi = \mathcal{P^{-1}_c}(\bf{m}) \tag{13}
Ψ=Pc−1(m)(13)
等式13表示从图像点
m
=
(
u
,
v
)
T
m = (u, v)^T
m=(u,v)T到入射光线方向
Ψ
=
(
θ
,
ψ
)
\Psi=(\theta, \psi)
Ψ=(θ,ψ)的入射。我们将
P
c
P_c
Pc写成一个复合的函数
P
c
=
A
∘
D
∘
F
\mathcal{P_c} = \mathcal{A} \circ \mathcal{D }\circ \mathcal{ F}
Pc=A∘D∘F,其中
F
\mathcal{F}
F表示等式7中从射线方向
Ψ
\Psi
Ψ到理想的卡氏坐标
x
=
(
x
,
y
)
T
\mathbf{x} = (x,y)^T
x=(x,y)T的变换。
D
\mathcal{D}
D表示从
x
\mathbf{x}
x到失真(畸变)坐标
x
d
=
(
x
d
,
y
d
)
T
\mathbf{x_d}=(x_d, y_d)^T
xd=(xd,yd)T的失真(畸变)映射,
A
\mathcal{A}
A表示等式11中的仿射变换。我们分解公式的原因是,对于变换
P
c
−
1
=
F
−
1
∘
D
−
1
∘
A
−
1
\mathcal{P}^{-1}_c = \mathcal{F}^{-1} \circ \mathcal{D}^{-1} \circ \mathcal{A}^{-1}
Pc−1=F−1∘D−1∘A−1,计算
F
−
1
\mathcal{F}^{-1}
F−1和
A
−
1
\mathcal{A}^{-1}
A−1非常直接,比较困难的部分是计算
D
−
1
\mathcal{D}^{-1}
D−1。
给一个点 x d \mathbf{x_d} xd,找 x = D − 1 ( x d ) \mathbf {x} = \mathcal{D}^{-1}(\mathbf{x_d}) x=D−1(xd)等价于计算 x = x d − s \mathbf{x = x_d -s} x=xd−s偏移量 s \mathbf{s} s,其中
s
=
S
(
Ψ
)
=
Δ
r
(
θ
,
ψ
)
u
r
(
ψ
)
+
Δ
t
(
θ
,
ψ
)
u
ψ
(
ψ
)
(14)
\mathbf{s} = \mathcal{S}(\Psi) = \Delta_r(\theta, \psi)\mathbf{u_r}(\psi) + \Delta_t(\theta, \psi)\mathbf{u_{\psi}}(\psi) \tag{14}
s=S(Ψ)=Δr(θ,ψ)ur(ψ)+Δt(θ,ψ)uψ(ψ)(14)
进一步说,我们可以写
S
(
Ψ
)
≡
(
S
∘
F
−
1
)
(
x
)
\mathcal{S(\Psi)}\equiv(\mathcal{S} \circ\mathcal{F}^{-1})(\mathbf{x})
S(Ψ)≡(S∘F−1)(x),并通过
S
∘
F
−
1
\mathcal{S} \circ\mathcal{F}^{-1}
S∘F−1的一阶泰勒公式在
x
d
\mathbf{x_d}
xd附近近似偏移量
s
\mathbf{s}
s得到
s
≃
(
S
∘
F
−
1
)
(
x
d
)
+
∂
(
S
∘
F
−
1
)
∂
x
(
x
d
)
(
x
−
x
d
)
=
S
(
Ψ
d
)
−
∂
S
∂
Ψ
(
∂
F
∂
Ψ
(
Ψ
d
)
)
−
1
s
\mathbf{s}\simeq(\mathcal{S} \circ\mathcal{F}^{-1})(\mathbf{x_d})+\frac{\partial (\mathcal{S} \circ\mathcal{F}^{-1}) }{\partial \mathbf{x}}(\mathbf{x_d})(\mathbf{x -x_d})= \mathcal{S}(\Psi_d) - \frac{\partial \mathcal{S}}{\partial \Psi}(\frac{\partial \mathcal{F}}{\partial \Psi}(\Psi_d))^{-1} \quad \mathbf{s}
s≃(S∘F−1)(xd)+∂x∂(S∘F−1)(xd)(x−xd)=S(Ψd)−∂Ψ∂S(∂Ψ∂F(Ψd))−1s
其中,
Ψ
d
=
F
−
1
(
x
d
)
\Psi_d = \mathcal{F}^{-1}(\mathbf{x_d})
Ψd=F−1(xd)可以被数值化估计,因此,我们可以从等式15计算偏移量
s
\mathbf{s}
s
s ≃ ( I + ∂ S ∂ Ψ ( Ψ d ) ( ∂ F ∂ Ψ ( Ψ d ) ) − 1 ) − 1 S ( Ψ d ) (15) \mathbf{s} \simeq(I + \frac{\partial \mathcal{S}}{\partial \Psi}(\Psi_d)(\frac{\partial \mathcal{F}}{\partial \Psi}(\Psi_d) )^{-1})^{-1} \mathcal{S}(\Psi_d) \tag{15} s≃(I+∂Ψ∂S(Ψd)(∂Ψ∂F(Ψd))−1)−1S(Ψd)(15)
其中,雅克比
∂
S
∂
Ψ
\frac{\partial \mathcal{S}}{\partial \Psi}
∂Ψ∂S和
∂
F
∂
Ψ
\frac{\partial \mathcal{F}}{\partial \Psi}
∂Ψ∂F从
等式14和等式7中计算,因此,最终
D
−
1
(
x
d
)
≃
x
d
−
(
I
+
(
∂
S
∂
Ψ
∘
F
−
1
)
(
x
d
)
(
(
∂
F
∂
Ψ
∘
F
−
1
)
(
x
d
)
)
−
1
)
−
1
(
S
∘
F
−
1
)
(
x
d
)
(16)
\mathcal{D}^{-1}(\mathbf{x_d}) \simeq\mathbf{x_d}-(I+(\frac{\partial \mathcal{S }}{\partial \Psi}\circ \mathcal{F}^{-1})(\mathbf{x_d})((\frac{\partial \mathcal{F}}{\partial \Psi}\circ\mathcal{F}^{-1})(\mathbf{x_d}))^{-1})^{-1} (\mathcal{S} \circ\mathcal{F}^{-1})(\mathbf{x_d}) \tag{16}
D−1(xd)≃xd−(I+(∂Ψ∂S∘F−1)(xd)((∂Ψ∂F∘F−1)(xd))−1)−1(S∘F−1)(xd)(16)
看起来非对称畸变(失真)函数
D
\mathcal{D}
D的一阶近似在实践中是成立的,因为回投模型的误差通常比正向模型的校准精度小几度,具体见第5节。
III. 投影模型矫正
传统的相机标定的方法以透视投影模型开始点,然后应用失真(畸变)项补偿。然而,对于鱼眼镜头这不是一个有效的方法,因为当 θ \theta θ接近 π / 2 \pi/2 π/2时,透视投影点到无线远处,这就导致了不可能用传统的畸变模型去除这种奇点。因此,我们将校准方法基于更通用的模型(即等式6)。
我们将等式6的投影模型与Mičušı́k16提出的针对鱼眼镜头的两个两参数模型作比较。
r
=
a
b
sin
(
b
θ
)
(
M
1
)
r
=
a
−
a
2
−
4
b
θ
2
2
b
θ
(
M
2
)
r = \frac{a}{b}\sin(b\theta) \quad (M1)\\ r = \frac{a-\sqrt{a^2-4b\theta^2}}{2b\theta} \quad (M2)
r=basin(bθ)(M1)r=2bθa−a2−4bθ2(M2)
图2中我们画了(1)–(5)投影曲线(即图一中的5个投影)以及模型
M
1
,
M
2
,
P
3
M1, M2, P3
M1,M2,P3 的最小二乘的近似,其中
P
3
P_3
P3指的是等式6的前两项。这里我们假设
f
=
200
f=200
f=200 pixels,这对于真实的相机也是一个合理的值。投影值在
[
0
,
θ
m
a
x
]
[0, \theta_{max}]
[0,θmax]之间近似,其中
θ
m
a
x
\theta_{max}
θmax分别取
6
0
∘
,
11
0
∘
,
11
0
∘
,
11
0
∘
,
9
0
∘
60^{\circ},110^{\circ},110^{\circ},110^{\circ},90^{\circ}
60∘,110∘,110∘,110∘,90∘。在
[
0
,
θ
m
a
x
]
[0,\theta_{max}]
[0,θmax]用0.1的步长离散化,
M
1
,
M
2
M1, M2
M1,M2用Levenberg-Marquardt拟合。从图2可以看出,
M
1
M1
M1对所有的透视投影和立体投影都不适合,
M
2
M2
M2对于正交投影不精确。
TABLE I中列列出了每个模型的最大误差,即图2中每个期望的曲线与近似曲线之间的最大水平距离。这里我们也模拟了
P
9
P9
P9,也就是多项式的前5项。可以看出模型
P
3
P3
P3在所有的2参数模型中有很好的性能,所有投影曲线的亚像素精度需要模型
P
9
P9
P9,这些结果表明,我们相机径向对称部分被很好的矫正了。
IV. 标定一般模型
接下来我们将要描述估计一般相机模型的步骤。标定方法基于观察已知位置中的控制点的平面物体(如棋盘格),与以前的方法相比,鱼眼镜头可能具有大于 18 0 ∘ 180^{\circ} 180∘的视野,也可以通过观察平面图案来标定。此外,如果使用圆形的控制点可以获得更高的精度,如IV-B中所述。
标定算法
标定程序包含下面所述的四个步骤。假设有N个视图的M个控制点,每一个视图有一个旋转矩阵
R
j
R_j
Rj和一个平移矩阵
t
j
t_j
tj描述相机相对于标定板平面的位置
X
c
=
R
j
X
+
t
j
,
j
=
1
,
.
.
.
,
N
(17)
\mathbf{X_c = R_jX+t_j}, \quad j=1,...,N \tag{17}
Xc=RjX+tj,j=1,...,N(17)
我们选择一个在X-Y平面内的标定平面,控制点
i
i
i的坐标用
X
i
=
(
X
i
,
Y
i
,
0
)
T
\mathbf{X}^i = (X^i, Y^i, 0)^T
Xi=(Xi,Yi,0)T表示。在标定平面上对应的齐次坐标表示为
x
p
i
=
(
X
i
,
Y
i
,
1
)
T
\mathbf{x}^i_p=(X^i, Y^i, 1)^T
xpi=(Xi,Yi,1)T,被
j
j
j视图观察的坐标表示为
m
j
i
=
(
u
j
i
,
v
j
i
)
\mathbf{m}^i_j=(u^i_j, v_j^i)
mji=(uji,vji)。标定程序的前三个步骤仅仅涉及到相机的6个内部参数
p
6
=
(
k
1
,
k
2
,
m
u
,
m
v
,
u
0
,
v
0
)
\mathbf{p}_6=(k_1,k_2, m_u, m_v,u_0,v_0)
p6=(k1,k2,mu,mv,u0,v0)。仅在最后一步插入完整模型的附加参数。
步骤一:内部参数的初始化
k 1 , k 2 k_1,k_2 k1,k2的初始化猜想可以通过拟合 r = k 1 θ + k 2 θ 3 r = k_1\theta+k_2\theta^3 r=k1θ+k2θ3拟合期望(2)–(4)的投影获得,其中制造商标称的焦距 f f f值和 θ m a x \theta_{max} θmax是已知的(可以根据手册查)。我们也可以获得传感器平面的最大图像半径 r m a x = k 1 θ m a x + k 2 θ 3 r_{max}=k_1\theta_{max}+k_2\theta^3 rmax=k1θmax+k2θ3。
带有圆形图像的鱼眼镜头,实际的图像仅仅填充图像帧内的圆形区域,在像素坐标中,该圆是个椭圆,它的参数可以被估计。
(
u
−
u
0
a
)
2
+
(
v
−
v
0
b
)
2
=
1
,
(\frac{u-u_0}{a})^2+(\frac{v-v_0}{b})^2 = 1,
(au−u0)2+(bv−v0)2=1,
因此,我们获取剩余参数
m
u
,
m
v
,
u
0
,
v
0
m_u,m_v, u_0,v_0
mu,mv,u0,v0的初始化猜想,其中
m
u
=
a
/
r
m
a
x
m_u =a/r_{max}
mu=a/rmax,
m
v
=
b
/
r
m
a
x
m_v=b/r_{max}
mv=b/rmax。对于full-frame相机镜头,最好的方法是将主点放在图像中心,并使用像素的维度(尺寸)值来获取
m
u
m_u
mu和
m
v
m_v
mv的初始值。
步骤二:反投影和单应性的计算
利用内部参数 p 6 = ( k 1 , k 2 , m u , m v , u 0 , v 0 ) T \mathbf{p6}=(k_1,k_2,m_u,m_v, u_0,v_0)^T p6=(k1,k2,mu,mv,u0,v0)T将观察点 m j i \mathbf{m}^{i}_j mji反投影到以相机源点为中心的单位球面上(参见图一(b)),球面上的点用 x ~ j i \tilde{\mathbf{x}}^i_j x~ji表示。由于标定平面上的点和单位球面上点的映射为中心投影,所以一个平面的单应性 H j \mathbf{H}_j Hj有 s x ~ j i = H j x p i s\tilde{\mathbf{x}}^i_j=\mathbf{H}_j\mathbf{x}^i_p sx~ji=Hjxpi。
对于每一个视图 j j j单应性 H j \mathbf{H}_j Hj通过下面计算:
(i) 首先通过计算归一化图像坐标坐标反投影控制点
(
x
j
i
y
j
i
)
=
[
1
/
m
u
0
0
1
/
m
v
]
(
u
j
i
−
u
0
v
j
i
−
v
0
)
\begin{pmatrix} x^i_j \\ y^i_j \end{pmatrix} = \begin{bmatrix} 1/m_u & 0 \\ 0 & 1/m_v \end{bmatrix} \begin{pmatrix} u^i_j-u_0 \\ v^i_j-v_0 \end{pmatrix}
(xjiyji)=[1/mu001/mv](uji−u0vji−v0)
变换成极坐标表示
(
r
j
i
,
ψ
j
i
)
=
^
(
x
j
i
,
y
j
i
)
(r^i_j,\psi^i_j)\hat{=}(x^i_j, y^i_j)
(rji,ψji)=^(xji,yji),最后解
θ
j
i
\theta^i_j
θji从立方等式
k
2
(
θ
j
i
)
+
k
1
θ
j
i
−
r
j
i
=
0
k_2(\theta^i_j)+k_1\theta^i_j-r^i_j=0
k2(θji)+k1θji−rji=0。
(ii) 设置
x
~
j
i
=
(
sin
ψ
j
i
sin
θ
j
i
,
cos
ψ
j
i
sin
θ
j
i
,
cos
θ
j
i
)
T
\tilde{\mathbf{x}}^i_j=(\sin\psi^i_j\sin\theta^i_j,\cos\psi^i_j\sin\theta^i_j,\cos\theta^i_j)^T
x~ji=(sinψjisinθji,cosψjisinθji,cosθji)T。
(iii)使用匹配点
x
^
j
i
⟺
x
p
i
\hat{\mathbf{x}}^i_j\iff\mathbf{x}^i_p
x^ji⟺xpi的数据归一化线性算法计算
H
j
\mathbf{H}_j
Hj的初始化估计。将
x
^
j
i
\hat{\mathbf{x}}^i_j
x^ji定义为在
H
j
\mathbf{H}_j
Hj条件下的精确图像,有
x
^
j
i
=
H
j
x
p
i
/
∣
∣
H
j
x
p
i
∣
∣
\hat{\mathbf{x}}^i_j=\mathbf{H}_j\mathbf{x}^i_p/||\mathbf{H}_j\mathbf{x}^i_p||
x^ji=Hjxpi/∣∣Hjxpi∣∣。
(iv) 最小化
∑
i
sin
2
α
j
i
\sum_i\sin^2\alpha^i_j
∑isin2αji精确单应性
H
j
\mathbf{H}_j
Hj,其中
α
j
i
\alpha^i_j
αji是
x
~
j
i
\tilde{\mathbf{x}}^i_j
x~ji和
x
^
j
i
\hat{\mathbf{x}}^i_j
x^ji之间的夹角。
步骤三:外部参数的初始化
相机外参的初始值从单应矩阵
H
j
\mathbf{H_j}
Hj提取。我们有
s
x
~
j
i
=
[
R
j
t
j
]
(
X
i
Y
i
0
1
)
=
[
r
j
1
r
j
2
t
j
]
(
X
i
Y
i
1
)
s\tilde{\mathbf{x}}^i_j =\mathbf{ [R_j \quad t_j]} \begin{pmatrix} X^i \\ Y^i \\ 0 \\ 1 \end{pmatrix}= \begin{bmatrix} \mathbf{r^1_j} & \mathbf{r^2_j} & \mathbf{t_j} \end{bmatrix} \begin{pmatrix} X^i \\ Y^i \\ 1 \end{pmatrix}
sx~ji=[Rjtj]⎝⎜⎜⎛XiYi01⎠⎟⎟⎞=[rj1rj2tj]⎝⎛XiYi1⎠⎞
上面的等式表明在差了一个尺度的情况下
H
j
=
[
r
j
1
r
j
2
t
j
]
\mathbf{H_j = [r^1_j \quad r^2_j \quad t_j]}
Hj=[rj1rj2tj],进一步有
r
j
1
=
λ
j
h
j
1
,
r
j
2
=
λ
j
h
j
2
,
r
j
3
=
r
j
1
×
r
j
2
,
t
j
=
λ
j
h
j
3
\mathbf{r^1_j=\lambda_jh^1_j,\qquad r^2_j=\lambda_jh^2_j,\qquad r^3_j=r^1_j \times r^2_j,\qquad t_j=\lambda_jh^3_j}
rj1=λjhj1,rj2=λjhj2,rj3=rj1×rj2,tj=λjhj3
其中
λ
j
=
s
i
g
n
(
H
j
3
,
3
)
/
∥
h
j
1
∥
\lambda_j = sign(H^{3,3}_j)/\Vert{\mathbf{h^1_j}}\Vert
λj=sign(Hj3,3)/∥hj1∥。由于估计误差的存在,获取的旋转矩阵不可能是正交的,因此,我们使奇异值分解计算在Frobenius 范数上最接近的正交矩阵并把它作为每个
R
j
\mathbf{R_j}
Rj的初始值猜想。
步骤四:最小化重投影误差
如果使用完整的模型(full model)
p
23
p_{23}
p23和
p
9
p_9
p9,则在此阶段将额外相机参数初始化为0。由于我们有内参和外参的估计,利用等式(17),(7)或(10),(11)计算每个相机的图像函数
P
j
\mathcal{P_j}
Pj,其中控制点被投影到
m
~
j
i
=
P
j
(
X
i
)
\mathbf{\tilde{m}^i_j}=\mathcal{P_j}(\mathbf{X^i})
m~ji=Pj(Xi)。相机参数通过测量值和模型的控制点之间的最小pingfang根和精确并使用Levenberg-Marquardt算法优化。
∑
j
=
1
N
∑
i
=
0
M
d
(
m
j
i
,
m
~
j
i
)
2
(18)
\sum_{j=1}^N\sum_{i=0}^Md(\mathbf{m^i_j,\tilde{m}^i_j})^2 \tag{18}
j=1∑Ni=0∑Md(mji,m~ji)2(18)
B. 圆形控制点的修改
为了实现精确校准,我们使用了带有白色圆圈黑色背景的校准平面,因为投影圆的质心可以达到亚像素级精度17。然而,在这个配置中投影圆的质心不是原始圆质心的图像。因此,由于式(18)中的
m
j
i
\mathbf{m^i_j}
mji是测量的质心,我们不应该把它作为
m
~
j
i
\mathbf{\tilde{m}^i_j}
m~ji质心。
为了避免上面的问题,我们提出了求解投影圆中心的方法,参数化以
(
X
0
,
Y
0
)
(X_0,Y_0)
(X0,Y0)为圆心,半径为
R
R
R的内部点
X
(
ρ
,
α
)
=
(
X
0
+
ρ
sin
α
,
Y
0
+
ρ
cos
α
,
0
)
T
\mathbf{X}(\rho,\alpha)=(X_0+\rho\sin\alpha, Y_0+\rho\cos\alpha, 0)^T
X(ρ,α)=(X0+ρsinα,Y0+ρcosα,0)T。给定相机参数,我们通过估计等式(19)得到圆的质心
m
~
\mathbf{\tilde{m}}
m~
m
~
=
∫
0
R
∫
0
2
π
m
~
(
ρ
,
α
)
∣
d
e
t
J
(
ρ
,
α
)
∣
d
α
d
ρ
∫
0
R
∫
0
2
π
∣
d
e
t
J
(
ρ
,
α
)
∣
d
α
d
ρ
\mathbf{\tilde{m}}=\frac{\int^R_0\int^{2\pi}_0\mathbf{\tilde{m}}(\rho,\alpha)|det\mathbf{J}(\rho,\alpha)|d\alpha d \rho}{\int^R_0\int^{2\pi}_0|det\mathbf{J}(\rho,\alpha)|d\alpha d \rho}
m~=∫0R∫02π∣detJ(ρ,α)∣dαdρ∫0R∫02πm~(ρ,α)∣detJ(ρ,α)∣dαdρ
其中 m ~ ( ρ , α ) = P ( X ( ρ , α ) ) \mathbf{\tilde{m}}(\rho,\alpha)=\mathcal{P}(\mathbf{X}(\rho,\alpha)) m~(ρ,α)=P(X(ρ,α))和 J ( ρ , α ) \mathbf{J}(\rho,\alpha) J(ρ,α)是组合函数 P ∘ X \mathcal{P}\circ\mathbf{X} P∘X的雅克比。雅克比的解析解是一项相当繁重的工作,但可以通过Maple等数学软件计算。
V. 标定实验
A. 传统和广角镜头相机
本文提出的方法与Heikki2中的使用的相机模型作对比,该模型具有4个畸变(失真)参数的零倾斜(skew-zero)的针孔模型,用 δ 8 \mathbf{\delta_8} δ8。
Brown, D. C.: “Close-range camera calibration”, Photogrammetric Engineering, 37(8):855-866, 1971. ↩︎
Heikkil¨ a, J.: “Geometric camera calibration using circular control points”, TPAMI, 22(10):1066-1077, 2000. ↩︎ ↩︎ ↩︎
Swaminathan, R. and Nayar, S. K.: “Nonmetric calibration of wide-angle lenses and polycameras”, TPAMI, 22(10):1172-1178, 2000. ↩︎
Miyamoto, K.: “Fish eye lens”, Journal of the Optical Society of America, 54(8):10601061, 1964. ↩︎
Bräuer-Burchardt, C. and Voss, K: “A new algorithm to correct fish-eye- and strong
wide-angle-lens-distortion from single images”, Proc. ICIP, pp. 225-228, 2001. ↩︎Devernay, F. and Faugeras, O.: “Straight lines have to be straight”, Machine Vision and
Applications, 13(1):14-24, 2001. ↩︎Basu, A. and Licardie, S.: “Alternative models for fish-eye lenses”, Pattern Recognition Letters, 16:433-441, 1995. ↩︎
Br¨auer-Burchardt, C. and Voss, K: “A new algorithm to correct fish-eye- and strong wide-angle-lens-distortion from single images”, Proc. ICIP, pp. 225-228, 2001 ↩︎
Devernay, F. and Faugeras, O.: “Straight lines have to be straight”, Machine Vision and Applications, 13(1):14-24, 2001. ↩︎
Shah, S. and Aggarwal, J.: “Intrinsic parameter calibration procedure for a (highdistortion) fish-eye lens camera with distortion model and accuracy estimation”, Pattern Recognition, 29(11):1775-1788, 1996. ↩︎
Bakstein, H. and Pajdla, T.: “Panoramic Mosaicing with a 180◦Field of View Lens”, Proc. IEEE Workshop on Omnidirectional Vision, pp. 60-67, 2002. ↩︎
Miˇcuˇ s´ık, B. and Pajdla, T.: “Estimation of omnidirectional camera model from epipolar geometry”, Proc. CVPR, pp. 485-490, 2003. ↩︎ ↩︎
Thirthala, S. and Pollefeys, M.: “Multi-view geometry of 1D radial cameras and its application to omnidirectional camera calibration”, Proc. ICCV , pp. 1539-1546, 2005. ↩︎ ↩︎
Claus, D. and Fitzgibbon, A. W.: “A rational function lens distortion model for general cameras”, Proc. CVPR, pp. 213-219, 2005. ↩︎ ↩︎
Barreto,J.P.andDaniilidis,K.:“Fundamentalmatrixforcameraswithradialdistortion”, Proc. ICCV ↩︎
Mičušı́k, B.: Two-View Geometry of Omnidirectional Cameras, PhD. Thesis, Czech
Technical University, 2004. ↩︎Heikkilä, J. and Silvén, O.: “Calibration procedure for short focal length off-the-shelf CCD cameras”, Proc. ICPR, pp. 166-170, 1996. ↩︎