A Flexible New Technique for Camera Calibration
我们提出了一个灵活、容易标定的新技术,它非常适合于在不需要特定已知的三维几何结构或者计算视觉中使用。该技术仅仅需要相机观察不同方向(至少两个方向)的平面图案(如棋盘格)。相机或者平面图案可以可以自由移动,运动的方式不需要知道。镜头的畸变使用径向畸变模型。提出的求解方式利用闭合解求解,然后以最大似然估计为基础的非线性优化方法提升参数精度。计算机模拟数据和真实的数据用来测试提出的标定技术,并取得了很好
目录
摘要
我们提出了一个灵活、容易标定的新技术,它非常适合于在不需要特定已知的三维几何结构或者计算视觉中使用。该技术仅仅需要相机观察不同方向(至少两个方向)的平面图案(如棋盘格)。相机或者平面图案可以自由移动,运动的方式不需要知道。镜头的畸变使用径向畸变模型。提出的求解方式利用闭合解求解,然后以最大似然估计为基础的非线性优化方法提升参数精度。计算机模拟数据和真实的数据用来测试提出的标定技术,并取得了很好的效果。与经典的技术,使用两个或三个正交的平面相比,本文提出的方法非常灵活,容易使用。它将3D计算机视觉从实验室环境一步推进到现实世界。
paper:A Flexible New Technique for Camera Calibration
author:Zhengyou Zhang
1 Motivations
在三维计算机视觉中,从2维图像中提取度量信息,相机标定是不可缺少的一步。从摄影测量界开始已做了大量工作1 2 3,最近在计算机视觉中更多4 5 6 7 8 9 10 11 。 我们可以将这些简单分为两类:摄影测量标定和自标定(photogrammetric calibration and self-calibration.)。
- 摄影测量标定:相机标定是通过观察已知的带有很高精度的三维空间几何结构的对象完成。该标定可以非常有效完成12 。 该标定对象通常由2个或3个非常精确的正交平面组成。有时,一个精确的平面的平移是已知的9 。 这些方法需要昂贵的校准设备和精心设置。
- 自标定: 此类别中的技术不使用任何校准对象,仅通过静态场景中相机的移动并利用它的图像信息就能够提供内部参数的2个约束13 14 。因此,如果图像由具有固定内部参数的相同相机拍摄,则三个图像之间的对应关系足以恢复内部和外部参数,这允许我们重建具有相似性的3-D结构15 16。虽然这种方法非常灵活,但尚未成熟17。 由于要估算的参数很多,我们无法始终获得可靠的结果。
存在其他技术:正交方向的消失点18 19和纯旋转校准20 21。
我们目前的研究主要集中在桌面视觉系统(DVS),因为使用DVS的潜力很大,相机变得便宜且无处不在。DVS针对的是普通大众,他们不是计算机视觉方面的专家。 典型的计算机用户将不时地执行视觉任务,因此不愿意为昂贵的设备投入资金。 因此,灵活性,稳健性和低成本是重要的。 本文中描述的相机标定技术是在考虑到这些因素的基础上开发的。
所提出的技术仅需要相机观察以几个(至少两个)不同方向的平面图案。该图案可以用激光打印机打印并附着在“感觉合理的”平面(例如,硬书封面)上。相机或者平面图案可以自由移动,运动的方式不需要知道。所提出的方法介于摄影测量标定和自标定之间,因为我们使用2D度量信息而不是3D或纯粹的隐式信息。计算机模拟数据和真实的数据用来测试提出的标定技术,并取得了很好的效果。与经典的技术相比,提出的技术相当的灵活,与自标定相比,它具有相当程度的稳健性。 我们相信这项新技术将3D计算机视觉从实验室环境向现实世界迈进了一步。
请注意,Bill Triggs 22最近从平面场景的至少5个视图开发了自标定技术,他的技术更加灵活,但是很难初始化。Liebowitz和Zisserman 19描述了一种使用度量信息(例如已知角度,两个相等但未知角度)和已知长度比率的平面透视图像的度量校正技术。 他们还提到,至少有三个这样的矫正平面可以校准内部摄像机参数,尽管没有显示实验结果。
本文的结构如下。 第2节描述了观察单个平面的基本约束。 第3节描述了标定程序。 我们从闭合解的开始,然后进行非线性优化。 还模拟了径向透镜畸变。 第4节研究了所提出的标定技术失败的配置。 在实践中很容易避免这种情况。 第5节提供了实验结果。 计算机模拟和实际数据都用于验证所提出的技术。 在附录中,我们提供了许多细节,包括估计模型平面与其图像之间的单应性的技术。
2 Basic Equations
我们通过观察单个平面来检查摄像机内部参数的约束。 我们从本文中使用的符号开始。
2.1 Notation
一个2D点用 m = [ u , v ] T \mathbf{m} = [u, v]^T m=[u,v]T表示,一个3D点用 M = [ X , Y , Z ] T \mathbf{M} = [X, Y, Z]^T M=[X,Y,Z]T表示。我们使用 x ^ \mathbf{\hat x} x^最后一个元素加1表示增广向量: m ~ = [ u , v , 1 ] T \mathbf{\widetilde{ m}} = [u, v, 1]^T m
=[u,v,1]T, M ~ = [ X , Y , Z , 1 ] \mathbf{\widetilde{ M}}=[X, Y, Z, 1] M
=[X,Y,Z,1]。相机通常使用针孔模型:一个3D点 M \mathbf{M} M和它的图像投影点 m \mathbf{m} m之间的关系是:
s m ~ = A [ R t ] M ~ (1) s \mathbf{\widetilde{ m} = A \begin{bmatrix} R & t \end{bmatrix} \widetilde{M} \tag{1} } sm
=A[Rt]M
(1)
其中s是任意的尺度因子, ( R , t ) (\mathbf{R, t}) (R,t)叫作外参,是世界坐标系统到相机坐标系统的旋转和平移关系, A \mathbf{A} A叫作相机的内参,写作
A = [ α γ u 0 0 β v 0 0 0 1 ] \mathbf{A} = \begin{bmatrix} \alpha & \gamma & u_0 \\ 0 & \beta & v_0 \\ 0 & 0 & 1 \end{bmatrix} A=⎣⎡α00γβ0u0v01⎦⎤
其中 u 0 , v 0 u_0, v_0 u0,v0是主点坐标, α , β \alpha, \beta α,β是图像u轴和v轴的尺度因子, γ \gamma γ是两个图像轴的偏斜参数。
我们使用缩写 A − T \mathbf{A}^{-T} A−T来表示 ( A − 1 ) T (\mathbf{A}^{-1})^T (A−1)T或 ( A T ) − 1 (\mathbf{A}^T)^{-1} (AT)−1。
2.2 Homography between the model plane and its image
为了不失一般性,我们假设模板平面在世界坐标系中表示为 Z = 0 Z=0 Z=0的平面。用 r i \mathbf{r}_i ri表示旋转矩阵 R \mathbf{R} R的第 i t h i^{th} ith列。从(1)中,我们有
s [ u v 1 ] = A = [ r 1 r 2 r 3 t ] [ X Y 0 1 ] = A [ r 1 r 2 t ] [ X Y 1 ] s\begin{bmatrix} u \\ v\\ 1 \end{bmatrix}= \mathbf{A} = \begin{bmatrix} \mathbf{r_1} & \mathbf{r_2} & \mathbf{r_3} & \mathbf{t} \end{bmatrix} \begin{bmatrix} X \\ Y \\ 0 \\ 1 \end{bmatrix} = \mathbf{A} \begin{bmatrix} \mathbf{r_1} & \mathbf{r_2} & \mathbf{t} \end{bmatrix} \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix} s⎣⎡uv1⎦⎤=A=[r1r2r3t]⎣⎢⎢⎡XY01⎦⎥⎥⎤=A[r1r2t]⎣⎡XY1⎦⎤
这里符号有些滥用,我们仍然使用 M M M表示模板平面上的一个点,但是 M = [ X , Y ] T \mathbf{M} = [X,Y]^T M=[X,Y]T,因为Z总是等于0。接着 M ~ = [ X , Y , 1 ] T \mathbf{\widetilde{M}}=[X,Y,1]^T M
=[X,Y,1]T。因此,一个模型点 M ~ \mathbf{\widetilde{M}} M
和它对应的 m \mathbf{m} m通过单应矩阵 H \mathbf{H} H关联:
s m ~ = H M ~ , w i t h H = A [ r 1 r 2 t ] (2) s \tilde{\mathbf{m}} = \mathbf{H}\tilde{\mathbf{M}}, \qquad with \qquad \mathbf{H} = \mathbf{A}\begin{bmatrix} \mathbf{r_1} & \mathbf{r_2} & \mathbf{t} \end{bmatrix} \tag{2} sm~=HM~,withH=A[r1r2t](2)
很清楚, H \mathbf{H} H是在尺度上相等的 3 × 3 3 \times 3 3×3矩阵。
2.3 Constraints on the intrinsic parameters
给定一个模型平面的图像,可以用来估计单应性矩阵(看附录A),使用 H = [ h 1 , h 2 , h 3 ] \mathbf{H=[h_1, h2, h3]} H=[h1,h2,h3]表示,从(2)中我们有:
[ h 1 h 2 h 3 ] = λ A [ r 1 r 2 t ] \begin{bmatrix} \mathbf{h1} & \mathbf{h2} & \mathbf{h3} \end{bmatrix} =\lambda \mathbf{A} \begin{bmatrix} \mathbf{r_1} & \mathbf{r_2} & \mathbf{t} \end{bmatrix} [h1h2h3]=λA[r1r2t]
其中, λ \lambda λ是任意的尺度因子,由于 r 1 , r 2 \mathbf{r_1, r_2} r1,r2是正交的,我们有
r 1 T r 2 = 0 h 1 T A − T A − 1 h 2 = 0 (3) \begin{aligned} \mathbf{r_1^T r_2} &= 0 \\ \mathbf{h^T_1A^{-T}A^{-1}h_2} &= 0 \tag{3} \end{aligned} r1Tr2h1TA−TA−1h2=0=0(3)
r 1 T r 1 = r 2 T r 2 h 1 T A − T A − 1 h 1 = h 2 T A − T A − 1 h 2 (4) \begin{aligned} \mathbf{r^T_1r_1 }&=\mathbf{ r^T_2r_2} \\ \mathbf{h^T_1A^{-T}A^{-1}h_1} &= \mathbf{h^T_2A^{-T}A^{-1}h_2} \tag{4} \end{aligned} r1Tr1h1TA−TA−1h1=r2Tr2=h2TA−TA−1h2(4)
给定一个单应性矩阵,对于内参矩阵来说有两个基本约束。因为一个单应矩阵有8个自由度,有六个是外参(3个旋转和3个平移),所以对于内部参数仅仅有2个约束。注意 A − T A − 1 \mathbf{A^{-T}A^{-1}} A−TA−1实际上描述了图像的绝对二次曲线16(absolute conic)。下一小节,我们从几何上解释。
2.4 Geometric Interpretation
(3),(4)跟绝对二次曲线相关,下面将推导如何得到(3),(4)。
依据惯例,很容易确定平面模型,在相机坐标系中使用下面平面方程描述:
[ r 3 r 3 T t ] T [ x y z ω ] = 0 \begin{bmatrix} \mathbf{r_3} \\ \mathbf{r^T_3t} \end{bmatrix}^T \begin{bmatrix} x \\ y \\ z \\ \omega \\ \end{bmatrix} = 0 [r3r3Tt]T⎣⎢⎢⎡xyzω⎦⎥⎥⎤=0
其中, ω = 0 \omega = 0 ω=0表示点在无限远处,否则 ω = 1 \omega = 1 ω=1。平面与平面在无限远处交于一条线,我们很容易看到 [ r 1 0 ] {\mathbf{r_1} \brack 0} [0r1]和 [ r 2 0 ] \mathbf{r_2} \brack 0 [0r2]是那条线上的特定点。因此在那条直线上的任何一个点可用这两个特定点的线性组合表示
x ∞ = a [ r 1 0 ] + b [ r 2 0 ] = [ a r 1 + b r 2 0 ] \mathbf{x_{\infin}} = a {\mathbf{r_1} \brack 0} + b {\mathbf{r_2} \brack 0} = {a \mathbf{r_1}+b \mathbf{r_2} \brack 0} x∞=a[0r1]+b[0r2]=[0ar1+br2]
现在,让我们计算上面直线与绝对二次曲线的交点。 根据定义,点 x ∞ \mathbf{x_{\infin}} x∞称为圆点,满足 x ∞ T x ∞ = 0 \mathbf{x^T_{\infin} x_{\infin}} = 0 x∞Tx∞=0,也就是
( a r 1 + b r 2 ) T ( a r 1 + b r 2 ) = 0 , 或 者 a 2 + b 2 = 0 (a\mathbf{r_1}+b\mathbf{r_2})^T(a\mathbf{r_1}+b\mathbf{r_2}) = 0, 或者 a^2+b^2=0 (ar1+br2)T(ar1+br2)=0,或者a2+b2=0
方程的解是 b = ± a i b = \pm ai b=±ai,其中 i 2 = − 1 i^2 =-1 i2=−1。也就是两个交点是
x ∞ = a [ r 1 + i r 2 0 ] \mathbf{x_{\infin}} = a {\mathbf{r_1}+i \mathbf{r_2} \brack 0} x∞=a[0r1+ir2]
然后得到他们在尺度上的图像平面上的投影为
m ~ ∞ = A ( r 1 ± i r 2 ) = h 1 ± i h 2 \mathbf{\widetilde{m}_{\infin}} = \mathbf{A(r_1 \pm} i\mathbf{r_2}) = \mathbf{h_1} \pm i \mathbf{h_2} m
∞=A(r1±ir2)=h1±ih2
点 m ~ \mathbf{\widetilde{m}} m
在绝对二次曲线的镜像上,由16 A − T A − 1 \mathbf{A^{-T}A^{-1}} A−TA−1描述,得到
( h 1 ± i h 2 ) T A − T A − 1 ( h 1 ± i h 2 ) = 0 (\mathbf{h_1 \pm}i \mathbf{h_2})^T\mathbf{A^{-T}A^{-1}}(\mathbf{h_1 \pm}i \mathbf{h_2}) = 0 (h1±ih2)TA−TA−1(h1±ih2)=0
因此,当要求实部和虚部均为零,得出(3)和(4)。
3 Solving Camera Calibration
这一节中详细描述如何有效的求解相机标定问题。我们以问题的解析解开始,接下来是以最大似然估计为基础的非线性最小二乘,最后我们将镜头的畸变考虑进来,给定解析解和非线性的解。
3.1 Closed-form solution
让
B = A − T A − 1 ≡ [ B 11 B 12 B 13 B 12 B 22 B 23 B 13 B 23 B 33 ] = [ 1 α 2 − γ α 2 β v 0 γ − u 0 β α 2 β − γ α 2 β γ 2 α 2 β 2 + 1 β 2 − γ ( v 0 γ − u 0 β ) α 2 β 2 − v 0 β 2 v 0 γ − u 0 β α 2 β − γ ( v 0 γ − u 0 β ) α 2 β 2 − v 0 β 2 ( v 0 γ − u 0 β ) 2 α 2 β 2 + v 0 2 β 2 + 1 ] (5) \mathbf{B = A^{-T}A^{-1}} \equiv \begin{bmatrix} B_{11} & B_{12} & B_{13} \\ B_{12} & B_{22} & B_{23} \\ B_{13} & B_{23} & B_{33} \end{bmatrix} \\ =\begin{bmatrix} \frac{1}{\alpha^2} & - \frac{\gamma}{\alpha^2 \beta} & \frac{v_0\gamma-u_0\beta}{\alpha^2\beta} \\ -\frac{\gamma}{\alpha^2\beta} & \frac{\gamma^2}{\alpha^2\beta^2}+\frac{1}{\beta^2} & -\frac{\gamma(v_0\gamma-u_0\beta)}{\alpha^2\beta^2}-\frac{v_0}{\beta^2} \\ \frac{v_0\gamma-u_0\beta}{\alpha^2\beta} & -\frac{\gamma(v_0\gamma-u_0\beta)}{\alpha^2\beta^2}-\frac{v_0}{\beta^2} & \frac{(v_0\gamma-u_0\beta)^2}{\alpha^2\beta^2}+\frac{v_0^2}{\beta^2}+1 \end{bmatrix} \tag{5} B=A−TA−1≡⎣⎡B11B12B13B12B22B23B13B23B33⎦⎤=⎣⎢⎡α21−α2βγα2βv0γ−u0β−α2βγα2β2γ2+β21−α2β2γ(v0γ−u0β)−β2v0α2βv0γ−u0β−α2β2γ(v0γ−u0β)−β2v0α2β2(v0γ−u0β)2+β2v02+1⎦⎥⎤(5)
注意 B \mathbf{B} B是对称的,定义6维向量 b \mathbf{b} b
b = [ B 11 , B 12 , B 22 , B 13 , B 23 , B 33 ] T (6) \mathbf{b} = [B_{11}, B_{12}, B_{22}, B_{13}, B_{23},B_{33}]^T \tag{6} b=[B11,B12,B22,B13,B23,B33]T(6)
设 H \mathbf{H} H的 i t h i^{th} ith列为 h i = [ h i 1 , h i 2 , h i 3 ] T \mathbf{h}_i=[h_{i1}, h_{i2},h_{i3}]^T hi=[hi1,hi2,hi3]T,然后有
h i T B h j = v i j T b (7) \mathbf{h}^T_i\mathbf{B}\mathbf{h}_j=\mathbf{v}^T_{ij}\mathbf{b} \tag{7} hiTBhj=vijTb(7)
其中
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 ] \mathbf{v}_{ij}=[h_{i1}h_{j1},h_{i1}h_{j2}+h_{i2}h_{j1},h_{i2}h_{j2}, \\ h_{i3}h_{j1}+h_{i1}h_{j3}, h_{i3}h_{j2}+h_{i2}h_{j3},h_{i3}h_{j3}] vij=[hi1hj1,hi1hj2+hi2hj1,hi2hj2,hi3hj1+hi1hj3,hi3hj2+hi2hj3,hi3hj3]
因此,结合(3),(4),(7)我们有
[ v 12 T ( v 11 − v 22 ) T ] b = 0 (8) \begin{bmatrix} \mathbf{v}^T_{12} \\ (\mathbf{v}_{11}-\mathbf{v}_{22})^T \end{bmatrix}\mathbf{b} = \mathbf{0} \tag{8} [v12T(v11−v22)T]b=0(8)
如果有n张观察模型平面(如棋盘格)的图像,通过堆叠n个(8)那样的等式,我们有
V b = 0 (9) \mathbf{Vb} = \mathbf{0} \tag{9} Vb=0(9)
其中 V \mathbf{V} V是 2 n × 6 2n\times6 2n×6的矩阵,如果n大于等于3,一般,我们将得到在尺度上的唯一解 b \mathbf{b} b。如果n=2,我们增加无偏约束(skewless constraint) γ = 0 \gamma=0 γ=0,即j将 [ 0 , 1 , 0 , 0 , 0 , 0 ] b = 0 [0, 1, 0, 0, 0, 0]\mathbf{b}=0 [0,1,0,0,0,0]b=0加到等式(9)中。(如果n=1,我们仅仅能解2个相机内参,例如 α , β \alpha, \beta α,β,同时假设 u 0 , v 0 u_0,v_0 u0,v0是已知的(例如在图像的中心), γ = 0 \gamma=0 γ=0,这确实是我们在中基于眼睛和嘴合理共面的事实确定头部姿势的方法23。),等式(9)的解是 V T V \mathbf{V^TV} VTV的最小奇异值对的奇异向量(等价地,与 V \mathbf{V} V的最小奇异值对应的右奇异向量)。
一旦估计了 b \mathbf{b} b,我们可以得到相机内参矩阵 A \mathbf{A} A,具体细节看附录 B \mathbf{B} B。一旦估计了 A \mathbf{A} A就可以轻松计算每个图像的外部参数(外参)。从(2)中,我们有
r 1 = λ A − 1 h 1 r 2 = λ A − 1 h 2 r 3 = r 1 × r 2 t = λ A − 1 h 3 \begin{alignedat}{4} \mathbf{r}_1 &= \lambda \mathbf{A}^{-1} \mathbf{h}_1 \\ \mathbf{r}_2 &= \lambda \mathbf{A}^{-1} \mathbf{h}_2 \\ \mathbf{r}_3 &= \mathbf{r}_1 \times \mathbf{r}_2 \\ \mathbf{t} &= \lambda \mathbf{A}^{-1} \mathbf{h}_3 \end{alignedat} r1r2r3t=λA−1h1=λA−1h2=r1×r2=λA−1h3
其中, λ = 1 / ∥ A − 1 h 1 ∥ = 1 / ∥ A − 1 h 2 ∥ \lambda=1/\Vert \mathbf{A}^{-1}\mathbf{h}_1 \Vert=1/\Vert \mathbf{A}^{-1}\mathbf{h}_2 \Vert λ=1/∥A−1h1∥=1/∥A−1h2∥。当然,由于数据中存在噪声,一般情况下所计算的矩阵 R = [ r 1 , r 2 , r 3 ] \mathbf{R}=[\mathbf{r}_1,\mathbf{r}_2,\mathbf{r}_3] R=[r1,r2,r3]不满足旋转矩阵的特性。附录 C \mathbf{C} C中描述了如何从一般的 3 × 3 3 \times 3 3×3矩阵估计最好的旋转矩阵的方法。
3.2 Maximum likelihood estimation
上面问题的解通过没有物理意义的最小代数距离来获得,下面 我们可以通过最大似然估计来精确它。
我们给定n张模型平面(例如棋盘格)的图像,每个图像有m个角点,假设这些图像点符合独立且均匀的噪声分布。最大似然估计可以通过下面的函数获得:
∑ i = 1 n ∑ j = 1 m ∥ m i j − m ^ ( A , R i , t i , M j ) ∥ 2 , (10) \sum_{i=1}^n \sum_{j=1}^m \Vert \mathbf{m}_{ij}-\hat{\mathbf{m}}(\mathbf{A,R}_i,\mathbf{t}_i,\mathbf{M}_j)\Vert^2, \tag{10} i=1∑nj=1∑m∥mij−m^(A,Ri,ti,Mj)∥2,(10)
其中 m ^ ( A , R i , t i , M j ) \hat{\mathbf{m}}(\mathbf{A,R}_i,\mathbf{t}_i,\mathbf{M}_j) m^(A,Ri,ti,Mj)根据等式(2)表示点 M j \mathbf{M}_j Mj在 i i i图像中的投影。旋转矩阵 R \mathbf{R} R被参数化为3个参数的向量 r \mathbf{r} r( r \mathbf{r} r平行于旋转轴,其长度大小等于旋转角度)组成, R , r \mathbf{R,r} R,r通过Rodrigues24公式转换。最小化(10)是一个非线性问题,可以通过Levenberg-Marquardt算法求解25,该方法需要初始化猜想 A \mathbf{A} A, R i , t i ∣ i = 1.. n {\mathbf{R}_i,\mathbf{t}_i|i=1..n} Ri,ti∣i=1..n,初始化猜想的值可以通过上一小节或得。
3.3 Dealing with radial distortion
直到现在我们都没有考虑相机的镜头畸变,但是,台式相机(desktop camera)通常会表现出明显的镜头畸变,尤其是径向畸变,在这一节中我们仅仅考虑径向畸变的前两项,更详细的模型,读者可以阅读26 1 3 11。根据报告文献3 927 ,畸变函数可能完全由径向分量主导,尤其是由第一项主导。也已经发现,任何更详细的建模不仅无济于事(与传感器量化相比可忽略不计),而且还将引起数值不稳定9 27。
让 ( u , v ) (u,v) (u,v)表示理想的(不可观测的、无畸变的)像素图像坐标, ( u ˘ , v ˘ ) (\breve u, \breve v) (u˘,v˘)表示对应的实际可观测到的图像坐标。理想点是根据针孔模型的模型点投影。类似的, ( x , y ) (x,y) (x,y)和 ( x ˘ , y ˘ ) (\breve x, \breve y) (x˘,y˘)分别是理想的(无畸变的)和实际的(带畸变的)归一化图像坐标。我们有3 27
x ˘ = x + x [ k 1 ( x 2 + y 2 ) + k 2 ( x 2 + y 2 ) 2 ] y ˘ = y + x [ k 1 ( x 2 + y 2 ) + k 2 ( x 2 + y 2 ) 2 ] \begin{alignedat}{2} \breve x = x + x[k_1(x^2+y^2)+k_2(x^2+y^2)^2] \\ \breve y= y + x[k_1(x^2+y^2)+k_2(x^2+y^2)^2] \end{alignedat} x˘=x+x[k1(x2+y2)+k2(x2+y2)2]y˘=y+x[k1(x2+y2)+k2(x2+y2)2]
其中 k 1 , k 2 k_1, k_2 k1,k2是畸变系数,径向畸变中心与主点相同。从 u ˘ = u 0 + α x ˘ + γ y ˘ \breve u = u_0+\alpha \breve x + \gamma \breve y u˘=u0+αx˘+γy˘, v ˘ = v 0 + β y ˘ \breve v = v_0 + \beta \breve y v˘=v0+βy˘以及假设 γ = 0 \gamma = 0 γ=0,我们有
u ˘ = u + ( u − u 0 ) [ k 1 ( x 2 + y 2 ) + k 2 ( x 2 + y 2 ) 2 ] (11) \breve u = u + (u-u_0)[k_1(x^2+y^2)+k_2(x^2+y^2)^2] \tag{11} u˘=u+(u−u0)[k1(x2+y2)+k2(x2+y2)2](11)
v ˘ = v + ( v − v 0 ) [ k 1 ( x 2 + y 2 ) + k 2 ( x 2 + y 2 ) 2 ] (12) \breve v= v + (v-v_0)[k_1(x^2+y^2)+k_2(x^2+y^2)^2] \tag{12} v˘=v+(v−v0)[k1(x2+y2)+k2(x2+y2)2](12)
迭代估计径向畸变 由于预期的径向畸变非常小,我们可以利用3.2节描述的技术估计其他的5个内部参数(通过简单的忽略径向畸变就可以做到,这样做也是合理的)。然后在估计相机内参(5个参数)后再去估计 k 1 , k 2 k_1, k_2 k1,k2,这样我们在估计 k 1 , k 2 k_1,k_2 k1,k2之前就会得到一个更加准确(理想)的坐标 ( u , v ) (u,v) (u,v),接着由(11)和(12),对于每个图像中的每个点我们有下面两个等式:
[ ( u − u 0 ) ( x 2 + y 2 ) ( u − u 0 ) ( x 2 + y 2 ) 2 ( v − v 0 ) ( x 2 + y 2 ) ( v − v 0 ) ( x 2 + y 2 ) 2 ] [ k 1 k 2 ] = [ u ˘ − u v ˘ − v ] \begin{bmatrix} (u-u_0)(x^2+y^2) & (u-u_0)(x^2+y^2)^2 \\ (v-v_0)(x^2+y^2) & (v-v_0)(x^2+y^2)^2 \end{bmatrix} \begin{bmatrix} k_1 \\ k_2 \end{bmatrix} =\begin{bmatrix} \breve u-u \\ \breve v - v \end{bmatrix} [(u−u0)(x2+y2)(v−v0)(x2+y2)(u−u0)(x2+y2)2(v−v0)(x2+y2)2][k1k2]=[u˘−uv˘−v]
给定 n n n个图像,每个图像 m m m个点,我们将等式堆叠在一起总共可得到 2 m n 2mn 2mn个等式,或者用矩阵形式可表示为 D k = d \mathbf{Dk = d} Dk=d,其中 k = [ k 1 , k 2 ] T \mathbf{k}=[k_1, k_2]^T k=[k1,k2]T。根据线性最小二乘得到:
k = ( D T D ) − 1 d (13) \mathbf{k=(D^TD)}^{-1}\mathbf{d} \tag{13} k=(DTD)−1d(13)
一旦估计了 k 1 , k 2 k_1, k_2 k1,k2,可利用(10)中的 m ^ ( A , R i , t i , M j ) \hat{\mathbf{m}}(\mathbf{A,R}_i,\mathbf{t}_i,\mathbf{M}_j) m^(A,Ri,ti,Mj)代替(11)和(12)进一步精确其他参数(相机内参,我们可以迭代这两个等式直至收敛(备注:即先求不带畸变的相机内参获得(u_0,v_0)->估计 k 1 , k 2 k_1,k_2 k1,k2->进一步优化相机内参得到(u_0,v_0)->进一步估计 k 1 , k 2 k_1,k_2 k1,k2->…,这样一直下去直至收敛)。
计算最大似然估计 实验中,我们发现上面的迭代技术很慢,一个很自然的扩展是最小化下面函数来求解参数集合:
∑ i = 1 n ∑ j = 1 m ∥ m i j − m ^ ( A , k 1 , k 2 , R i , t i , M j ) (14) \sum_{i=1}^n \sum_{j=1}^m \Vert \mathbf{m}_{ij}-\hat{\mathbf{m}}(\mathbf{A},k_1, k_2,\mathbf{R}_i,\mathbf{t}_i,\mathbf{M}_j) \tag{14} i=1∑nj=1∑m∥mij−m^(A,k1,k2,Ri,ti,Mj)(14)
根据(2), m ^ ( A , k 1 , k 2 , R i , t i , M j ) \hat{\mathbf{m}}(\mathbf{A},k_1, k_2,\mathbf{R}_i,\mathbf{t}_i,\mathbf{M}_j) m^(A,k1,k2,Ri,ti,Mj)表示第 i i i个图像的第 M j M_j Mj在遵循(11)和(12)畸变条件下的投影。(14)是非线性问题,可以通过Minpack28中的Levenberg-Marquardt算法求解。 R \mathbf{R} R参数化为三维向量 r \mathbf{r} r(详见3.2节), A \mathbf{A} A和 { R i , t i ∣ i = 1.. n } \{\mathbf{R}_i, \mathbf{t}_i|i=1..n\} {Ri,ti∣i=1..n}可以通过3.1节和3.2节描述的技术求解。 k 1 , k 2 k_1,k_2 k1,k2可以通过上一部分描述的方法获得,也可以简单的设置它们为0。
3.4 Summary
推荐的标定程序如下:
- 打印一个棋盘格并固定在平面(面板)上
- 在不同方向上移动面板(或者相机)采取模板(棋盘格)的一组图像
- 在图像中检测特征点(棋盘格点)
- 使用闭合解(3.1节描述)估计5个内部参数和所有的外部参数
- 根据(13)描述的线性最小二乘估计径向系数
- 最小化(14)精确所有的参数
4. Degenerate Configurations
本节我们研究估计相机内参时,增加图像的数量并不能增加更多的约束的配置。因为(3)和(4)的推导来源于旋转矩阵,如果 R 2 \mathbf{R}_2 R2和 R 1 \mathbf{R}_1 R1不相关,则图像2不能提供额外的约束。实际上,如果一个平面暗含纯平移,则 R 2 = R 1 \mathbf{R}_2 = \mathbf{R}_1 R2=R1,图像2对于相机的标定毫无帮助。接下来我们考虑更加复杂的配置。
命题 1. 如果模型平面(如棋盘格平面)在第二位置平行于其第一位置,则第二单应性不提供其他约束。
证明:一般情况下, R 2 \mathbf{R}_2 R2和 R 1 \mathbf{R}_1 R1通过绕z轴旋转相关,也就是
R 1 [ cos θ − sin θ 0 sin θ cos θ 0 0 0 1 ] = R 2 \mathbf{R}_1 \begin{bmatrix} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{bmatrix} = \mathbf{R}_2 R1⎣⎡cosθsinθ0−sinθcosθ0001⎦⎤=R2
其中 θ \theta θ是相关的旋转角,我们使用上标(1),(2)分别表示与图像1,2相关的向量。为了清晰我们有
h 1 ( 2 ) = λ ( 2 ) ( A r ( 1 ) cos θ + A r ( 2 ) sin θ ) = λ ( 2 ) λ ( 1 ) ( h 1 ( 1 ) cos θ + h 2 ( 1 ) sin θ ) h 2 ( 2 ) = λ ( 2 ) ( − A r ( 1 ) sin θ + A r ( 2 ) cos θ ) = λ ( 2 ) λ ( 1 ) ( − h 1 ( 1 ) sin θ + h 2 ( 1 ) cos θ ) \begin{alignedat}{2} \mathbf{h}^{(2)}_1 = \lambda^{(2)}(\mathbf{Ar}^{(1)} \cos \theta + \mathbf{Ar}^{(2)} \sin \theta)=\frac{\lambda^{(2)}}{\lambda^{(1)}}(\mathbf{h}^{(1)}_1 \cos \theta + \mathbf{h}^{(1)}_2 \sin \theta) \\ \mathbf{h}^{(2)}_2 = \lambda^{(2)}(-\mathbf{Ar}^{(1)} \sin \theta + \mathbf{Ar}^{(2)} \cos \theta)=\frac{\lambda^{(2)}}{\lambda^{(1)}}(-\mathbf{h}^{(1)}_1 \sin \theta + \mathbf{h}^{(1)}_2 \cos \theta) \end{alignedat} h1(2)=λ(2)(Ar(1)cosθ+Ar(2)sinθ)=λ(1)λ(2)(h1(1)cosθ+h2(1)sinθ)h2(2)=λ(2)(−Ar(1)sinθ+Ar(2)cosθ)=λ(1)λ(2)(−h1(1)sinθ+h2(1)cosθ)
接着图像2的约束(3)变成:
h 1 ( 2 ) T A − T A − 1 h 2 ( 2 ) = λ ( 2 ) λ ( 1 ) [ ( cos 2 θ − sin 2 θ ) ( h 1 ( 1 ) T A − T A − 1 h 2 ( 1 ) ) − cos θ sin θ ( h 1 ( 1 ) T A − T A − 1 h 1 ( 1 ) ) − h 2 ( 1 ) T A − T A − 1 h 2 ( 1 ) ] \mathbf{h}^{(2)^T}_1 \mathbf{A}^{-T} \mathbf{A}^{-1} \mathbf{h}^{(2)}_2 = \frac{\lambda^{(2)}}{\lambda^{(1)}}[(\cos^2 \theta - \sin^2 \theta)(\mathbf{h}^{(1)^T}_1 \mathbf{A}^{-T} \mathbf{A}^{-1} \mathbf{h}^{(1)}_2) \\ -\cos \theta \sin \theta(\mathbf{h}^{(1)^T}_1 \mathbf{A}^{-T} \mathbf{A}^{-1} \mathbf{h}^{(1)}_1)- \mathbf{h}^{(1)^T}_2 \mathbf{A}^{-T} \mathbf{A}^{-1} \mathbf{h}^{(1)}_2 ] h1(2)TA−TA−1h2(2)=λ(1)λ(2)[(cos2θ−sin2θ)(h1(1)TA−TA−1h2(1))−cosθsinθ(h1(1)TA−TA−1h1(1))−h2(1)TA−TA−1h2(1)]
这是 H 1 \mathbf{H}_1 H1提供的两个约束的线性组合。 同样,我们可以证明图像2的第二个约束也是 H 1 \mathbf{H}_1 H1提供的两个约束的线性组合,因此,我们没有从 H 2 \mathbf{H}_2 H2获得任何约束。
结果是不言而喻的,因为平行平面与平面在相同的圆形点处在无限远处相交,因此符合2.4节它们提供相同的约束。
实际上,避免简并配置非常容易:我们只需要将模型平面(如棋盘格)的方向从一个方向更改为另一个方向即可。
尽管如果模型平面进行纯平移,提出的技术将不起作用,但如果知道平移,则仍然可以进行相机校准。 请参考附录D。
附录D Camera Calibration Under Known Pure Translation
正如第4节所说的那样,如果模型平面仅做平移,那么提出的技术不起作用。然而,如果平移已知,相机标定是可能的,例如Tsai’s9的技术。从(2)有 t = α A − 1 h 3 \mathbf{t}=\alpha \mathbf{A}^{-1}\mathbf{h}_3 t=αA−1h3,其中 α = 1 A − 1 h 1 \alpha = \frac{1}{\mathbf{A}^{-1}\mathbf{h}_1} α=A−1h11,在 i , j i,j i,j两个位置间的平移可表示为
t i j = t i − t j = A − 1 ( α i h 3 i − α j h 3 j ) \mathbf{t}^{ij} = \mathbf{t}^{i} - \mathbf{t}^{j} = \mathbf{A}^{-1}(\alpha^{i}\mathbf{h}^{i}_3 - \alpha^{j}\mathbf{h}^{j}_3) tij=ti−tj=A−1(αih3i−αjh3j)
(请注意,尽管 H ( i ) \mathbf{H}^{(i)} H(i)和 H ( j ) \mathbf{H}^{(j)} H(j)均按其自身的比例因子估算,但可以使用纯平移的因子将它们重新按比例缩放为单个公共比例因子)。如果仅仅知道平移的方向,我们可以在 A \mathbf{A} A上得到两个约束。如果我们还知道平移的大小(magnitude),那我们又增加了一个约束在 A \mathbf{A} A上,那么两个平面就可以完成所有标定。
-
D. C. Brown. Clos-range camera calibration. Photogrammetric Engineering, 37(8):855–866,1971. ↩︎ ↩︎
-
B. Caprile and V. Torre. Using Vanishing Points for Camera Calibration. The International Journal of Computer Vision, 4(2):127–140, Mar. 1990. ↩︎
-
W. Faig. Calibration of close-range photogrammetry systems: Mathematical formulation. Pho-
togrammetric Engineering and Remote Sensing, 41(12):1479–1486, 1975. ↩︎ ↩︎ ↩︎ ↩︎ -
O. Faugeras, T. Luong, and S. Maybank. Camera self-calibration: theory and experiments. In
G. Sandini, editor, Proc 2nd ECCV, volume 588 of Lecture Notes in Computer Science, pages
321–334, Santa Margherita Ligure, Italy, May 1992. Springer-Verlag. ↩︎ -
O. Faugeras and G. Toscani. The calibration problem for stereo. In Proceedings of the IEEE
Conference on Computer Vision and Pattern Recognition, pages 15–20, Miami Beach, FL, June,1986. IEEE. ↩︎ -
S. Ganapathy. Decomposition of transformation matrices for robot vision. Pattern Recognition
Letters, 2:401–412, Dec. 1984. ↩︎ -
D. Gennery. Stereo-camera calibration. In Proceedings of the 10th Image Understanding Work-
shop, pages 101–108, 1979. ↩︎ -
S. J. Maybank and O. D. Faugeras. A theory of self-calibration of a moving camera. The
International Journal of Computer Vision, 8(2):123–152, Aug. 1992. ↩︎ -
R. Y. Tsai. A versatile camera calibration technique for high-accuracy 3D machine vision
metrology using off-the-shelf tv cameras and lenses. IEEE Journal of Robotics and Automa-
tion, 3(4):323–344, Aug. 1987. ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ -
G. Wei and S. Ma. A complete two-plane camera calibration method and experimental compar-
isons. In Proc. Fourth International Conference on Computer Vision, pages 439–446, Berlin,
May 1993. ↩︎ -
J. Weng, P. Cohen, and M. Herniou. Camera calibration with distortion models and accuracy
evaluation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(10):965–980,
Oct. 1992. ↩︎ ↩︎ -
O. Faugeras. Three-Dimensional Computer Vision: a Geometric Viewpoint. MIT Press, 1993. ↩︎
-
Q.-T. Luong. Matrice Fondamentale et Calibration Visuelle sur l’Environnement-Vers une plus
grande autonomie des systèmes robotiques. PhD thesis, Université de Paris-Sud, Centre d’Orsay,
Dec. 1992. ↩︎ -
S. J. Maybank and O. D. Faugeras. A theory of self-calibration of a moving camera. The
International Journal of Computer Vision, 8(2):123–152, Aug. 1992. ↩︎ -
R. I. Hartley. An algorithm for self calibration from several views. In Proceedings of the IEEE
Conference on Computer Vision and Pattern Recognition, pages 908–912, Seattle, WA, June1994. IEEE. ↩︎ -
Q.-T. Luong and O. Faugeras. Self-calibration of a moving camera from point correspondences
and fundamental matrices. The International Journal of Computer Vision, 22(3):261–289, 1997. ↩︎ ↩︎ ↩︎ -
S. Bougnoux. From projective to euclidean space under any practical situation, a criticism of
self-calibration. In Proceedings of the 6th International Conference on Computer Vision, pages
790–796, Jan. 1998. ↩︎ -
B. Caprile and V. Torre. Using Vanishing Points for Camera Calibration. The International
Journal of Computer Vision, 4(2):127–140, Mar. 1990. ↩︎ -
D. Liebowitz and A. Zisserman. Metric rectification for perspective images of planes. In Pro-
ceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 482–488,
Santa Barbara, California, June 1998. IEEE Computer Society. ↩︎ ↩︎ -
R. Hartley. Self-calibration from multiple views with a rotating camera. In J.-O. Eklundh, editor,
Proceedings of the 3rd European Conference on Computer Vision, volume 800-801 of Lecture
Notes in Computer Science, pages 471–478, Stockholm, Sweden, May 1994. Springer-Verlag. ↩︎ -
G. Stein. Accurate internal camera calibration using rotation, with analysis of sources of er-
ror. In Proc. Fifth International Conference on Computer Vision, pages 230–236, Cambridge,
Massachusetts, June 1995. ↩︎ -
B. Triggs. Autocalibration from planar scenes. In Proceedings of the 5th European Conference
on Computer Vision, pages 89–105, Freiburg, Germany, June 1998. ↩︎ -
I. Shimizu, Z. Zhang, S. Akamatsu, and K. Deguchi. Head pose determination from one image
using a generic model. In Proceedings of the IEEE Third International Conference on Automatic
Face and Gesture Recognition, pages 100–105, Nara, Japan, Apr. 1998. ↩︎ -
Faugeras. Three-Dimensional Computer Vision: a Geometric Viewpoint. MIT Press, 1993. ↩︎
-
O. Faugeras. Three-Dimensional Computer Vision: a Geometric Viewpoint. MIT Press, 1993. ↩︎
-
C. Slama, editor. Manual of Photogrammetry. American Society of Photogrammetry, fourth
edition, 1980. ↩︎ -
G. Wei and S. Ma. Implicit and explicit camera calibration: Theory and experiments. IEEE
Transactions on Pattern Analysis and Machine Intelligence, 16(5):469–480, 1994. ↩︎ ↩︎ ↩︎ -
J. More. The levenberg-marquardt algorithm, implementation and theory. In G. A. Watson,
editor, Numerical Analysis, Lecture Notes in Mathematics 630. Springer-Verlag, 1977. ↩︎
更多推荐
所有评论(0)