Camera Calibration 相机标定：原理简介（五）

来源：互联网发布：python 打包 linux 编辑：程序博客网时间：2024/05/07 06:41

5 基于2D标定物的标定方法

基于2D标定物的标定方法，原理与基于3D标定物相同，只是通过相机对一个平面进行成像，就可得到相机的标定参数，由于标定物为平面，本身所具有的约束条机，相对后者标定更为简单。经典算法为Z. Zhang(PAMI, 2000) A Flexible New Technique for Camera Calibration。其算法已经被收入Opencv(2004)，最常用的标定图案是棋盘格图案，如下图：

ChessBoard

5.1 单应性矩阵

对于2D标定平面，抑或称为标定板，不妨假设，平面上点的增广齐次向量[X,Y,Z,1]满足Z=0，同时对于旋转矩阵R的每一列表示为一个列向量记为ri，则根据相机的投影方程：

s ⎡ ⎣ ⎢ u v 1 ⎤ ⎦ ⎥ = A [r 1 r 2 r 3 t] ⎡ ⎣ ⎢ ⎢ ⎢ X Y 01 ⎤ ⎦ ⎥ ⎥ ⎥ = A [r 1 r 2 t] ⎡ ⎣ ⎢ X Y 1 ⎤ ⎦ ⎥ (1)

因为点[X,Y]T仍表示的是三维空间点坐标，只是由于标定使用平面的特殊性，将Z=0省略，因此我们仍然使用统一的M符号表示，与其对应M~=[X,Y,1]T表示该点的其次向量，则公式(1)可简记为：

s m ~ = H M ~ with H = A [r 1 r 2 t] (2)

H被称为单应性矩阵（Homography matrix）。

5.2 内参约束条件

对于单应性矩阵H，我们也将其按照列向量的方式表示：H=[h1 h2 h3]，则有：

[h 1 h 2 h 3] = λ A [r 1 r 2 t] (3)

其中，λ是一个任意的系数，因为r1,r2是正交向量，因此有：

r T 1 r 2 = h T 1 A - T A - 1 h 2 = 0 (4)

r T 1 r 1 = h T 1 A - T A - 1 h 1 = r T 2 r 2 = h T 2 A - T A - 1 h 2 (5)

对于一个单应性矩阵，公式(4,5)是两个内参的基本约束。单应性矩阵有9个元素，但可以由8个独立不相关的元素表示，也就是说投影变换有8个自由度，而我们只有6个外参元素（3个旋转元素和3个平移元素），因此两个内参约束条件的意义就在于此。在原理简介（三）中，我们已经了解到A−TA−1描述的就是IAC（Image of the absolute conic），后文将对此进行解释。

5.3 几何解释

现在让我们来分析公式(4,5)与绝对圆锥曲线的关系。在像空间坐标系中，存在下面的等式：

[r 3 r T 3 t] T ⎡ ⎣ ⎢ ⎢ ⎢ x y z w ⎤ ⎦ ⎥ ⎥ ⎥ = 0 (6)

当w=0时就是我们前面解释的改点位于无穷远处。我们想象标定板所在平面与该无穷远处平面相交于一点，则点[r10]和[r20]是相交直线上的两个特殊点，线上的其他点都可以用这两个点线性表示：

x \infty = a [r 1 0] + b [r 2 0] = [a r 1 + b r 2 0] (7)

现在让我们看一下上述的交线和绝对圆锥曲线的交点，原理简介（三）已经介绍点x∞满足：xT∞x∞=0，也即(ar1+br2)T(ar1+br2)=0 ⇒ a2+b2=0（r1与r2正交）。因此b=±ai, 其中i2=−1，公式(7)可写为：

x \infty = a [r 1 \pm i r 2 0] (8)

其在图像中的投影点为：

m \infty = A (r 1 \pm i r 2) = h 1 \pm i h 2 (9)

因为点m∞位于IAC上，有：

(h 1 \pm i h 2) T A - T A - 1 (h 1 \pm i h 2) = 0 (10)

因此等式(10)左边的实部与虚部都为0，也就是公式(4,5)这两条约束条件。

5.4 闭合解

现在开始讲解如何高效求解本方法的相机标定问题。做法是，首先获得分析解，然后初始估计值基于最大似然准则进行非线性优化，这些都将逐步进行讲解。

令：

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 2 0 β 2 + 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ (11)

因此可以看出，B是一个对称矩阵，由六个元素组成：

b = [B 11, B 12, B 13, B 22, B 23, B 33] T (12)

对于H=[h1,h2,h3]的第i列记为：hi=[hi1,hi2,hi3]T，则有：

h T i B h j = v T i j b (13)

其中，vij=[hi1hj1,hi1hj2+hi2hj1,hi2hj2,hi3hj1+hi1hj3,hi3hj2+hi2hj3,hi3hj3]T。因此，对于公式(4,5)可以表达为齐次方程：

[v T 12 (v 11 - v 22) T] b = 0 (14)

如果采集了n张图像，可以列出方程组：

V b = 0 (15)

其中，V是一个2n×6的矩阵。如果n⩾3，一般就能获得b的唯一解；如果n=2，我们可以利用偏斜度γ=0的约束条件，将[0,1,0,0,0,0]b=0添加到方程组(15)中；如果n=1那么就只能获得相机的内参。

一旦b确认后，就可以获得相机内参。将公式(11)中的矩阵B添加一个任意系数：B=λA−TA，可得到：

v 0 = (B 12 B 13 - B 11 B 23) / (B 11 B 22 - B 212) (16)

λ = B 33 - [B 213 + v 0 (B 12 B 13 - B 11 B 23)] / B 11 (17)

α = λ / B 11 - - - - - \sqrt (18)

β = λ B 11 / (B 11 B 12 - B 212) - - - - - - - - - - - - - - - - - \sqrt (19)

γ = - B 12 α 2 β / λ (20)

u 0 = γ v 0 / α - B 13 α 2 / λ (21)

内参矩阵获得后，根据公式(2)外参也就很快就可以得到：

r 1 = λ' A - 1 h 1, r 2 = λ' A - 1 h 2, r 3 = r 1 \times r 2, t = λ' A - 1 h 3 (22)

其中，λ′=1/∥A−1h1∥=1/∥A−1h2∥。

由于噪声的存在，这样求解的矩阵R一般并不具备旋转矩阵的特性，更好的求解方法是通过奇异值分解的方法，可以参考Z. Zhang的论文。

5.5 最大似然优化

与基于3D标定物的优化方法类似，这里仍然使用最大似然法进行优化，也就是基于噪声是不相关且独立分布的的假设，对于n张图像，其中包含m个点，优化方程为：

\sum i = 0 n \sum j = 0 m ∥ m i j - m^(A, R i, t i, M j) ∥ 2 (23)

其中，m^(A,Ri,ti,Mj)是三维点Mj在图像i中的投影点。非线性优化过程是，首先需要获得内参矩阵A的估算值，然后利用上面的描述的求解方法，获得{Ri,ti|i=1,…,n}，利用LM迭代优化方法，使得公式(23)的求和最小化。

5.6 镜头畸变

在原理简介（四）4.6中已经对镜头畸变模型进行了讲解，这里与之原理一样。以径向畸变为例，为了获得较为理想准确的像点坐标值，则有方程：

[(u - u 0) (x 2 + y 2) (v - v 0) (x 2 + y 2) (u - u 0) (x 2 + y 2) 2 (v - v 0) (x 2 + y 2) 2)] [k 1 k 2] = [u^- u v^- v] (24)

同样当m个点在n张图像中时，可以组成方程组Dk=d，其中k=[k1k2]，其线性最小二乘解为：

k = (D T D) - 1 D T d (25)

当获得(k1,k2)后，我们就可以将公式(23)的非线性优化调整为：

\sum i = 0 n \sum j = 0 m ∥ m i j - m^(A, k 1, k 2, R i, t i, M j) ∥ (26)

结束语：
基于2D标定物的标定方法流程为：

1 打印出一张标定图并贴到一个平面上；
2 通过移动相机或者标定平面采集不同位置、不同方向的标定板图像；
3 特征点检测；
4 估算内参（公式(16−21)），然后通过闭合解得到外参（公式(22)）；
5 通过最小二乘线性估算畸变系数；
6 使公式(26)的代数和最小，优化所有参数。

1 0