>>>文章索引<<<

相机标定（一）-原理及内参、外参

相机标定（二）-畸变校正，张正友标定法

相机标定（三）-相机成像模型

1 一些基本的方程推导

1.1 预定义

定义2D点为 m = [ u , v ] T m=[u,v]^T m=[u,v]T，3D点为 M = [ X , Y , Z ] T M=[X,Y,Z]^T M=[X,Y,Z]T。用上标~用来表示增广向量（在最后添加元素1，齐次形式）：

m ~ = [ u , v , 1 ] T \tilde{m}=[u,v,1]^T m~=[u,v,1]T

M ~ = [ X , Y , Z ] T \tilde{M}=[X,Y,Z]^T M~=[X,Y,Z]T

摄像机通常建模为针孔模型，则三维点 M M M及其在图像上的投影 m m m之间的关系可以表示为：

s m ~ = A [ R t ] M ~ ,其中 A = [ α c u 0 0 β v 0 0 0 1 ] (1) s\tilde{m}=A \begin{bmatrix} R & t \end{bmatrix} \tilde{M} \text{, 其中} A= \begin{bmatrix} \alpha & c & u_0 \\ 0 & \beta & v_0 \\ 0 & 0 & 1 \end{bmatrix} \tag{1} sm~=A[R​t​]M~,其中A=⎣⎡​α00​cβ0​u0​v0​1​⎦⎤​(1)

其中 s s s为一个任意的尺度因子（与拍摄的标定图像相关）； ( R , t ) (R,t) (R,t)为相机的外参矩阵，表示相机坐标系与世界坐标系之间的转换； A A A为相机的内参矩阵，其中 ( u 0 , v 0 ) (u_0,v_0) (u0​,v0​)为相机主点的坐标， α \alpha α和 β \beta β为图像 u u u和 v v v坐标轴上的尺度因子， c c c描述了图像两个坐标轴之间的偏度。（关于相机内、外参矩阵的推导可参考：链接）

本文中使用 A − T A^{-T} A−T来表示 ( A − 1 ) T (A^{-1})^T (A−1)T和 ( A T ) − 1 (A^T)^{-1} (AT)−1。

1.2 棋盘格模型与其图像之间的单应性矩阵

文中把世界坐标定义在棋盘格模型上，即模型平面与世界坐标系中 Z = 0 Z=0 Z=0的平面重合。我们定义旋转矩阵 R R R的第 i i i列为 r i r_i ri​，则由公式 ( 1 ) (1) (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 ] \begin{aligned} s \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} &= A \begin{bmatrix} r_1 & r_2 & r_3 & t \end{bmatrix} \begin{bmatrix} X \\ Y \\ 0 \\ 1 \end{bmatrix} \\ &= A \begin{bmatrix} r_1 & r_2 & t \end{bmatrix} \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix} \end{aligned} s⎣⎡​uv1​⎦⎤​​=A[r1​​r2​​r3​​t​]⎣⎢⎢⎡​XY01​⎦⎥⎥⎤​=A[r1​​r2​​t​]⎣⎡​XY1​⎦⎤​​

由于我们定义 Z = 0 Z=0 Z=0，所以棋盘格模型上的点为 M = [ X , Y ] T M=[X,Y]^T M=[X,Y]T，其齐次形式为 M ~ = [ X , Y , 1 ] T \tilde{M}=[X,Y,1]^T M~=[X,Y,1]T。这样，模型上的点 M M M与其图像上对应成像的点 m m m就可以通过一个单应性矩阵 H H H来联系：

s m ~ = H M ~ ,其中 H = A [ r 1 r 2 t ] (2) s\tilde{m}=H\tilde{M} \text{, 其中} H=A \begin{bmatrix} r_1 & r_2 & t \end{bmatrix} \tag{2} sm~=HM~,其中H=A[r1​​r2​​t​](2)

显然，这个 3 × 3 3\times3 3×3的单应性矩阵 H H H的定义依赖于尺度因子 s s s。

1.3 内参的约束条件

给定一幅棋盘格模型的图像，可以估计一个单应性矩阵。文中作者使用一种极大似然估计的方法（可参考链接）来求解。

假设 M i M_i Mi​为模型点， m i m_i mi​为图像点，则理想情况下它们满足公式 ( 2 ) (2) (2)，但实际上由于图像中存在噪声，无法满足公式 ( 2 ) (2) (2)。作者假设图像中的噪声满足均值为 0 0 0、协方差矩阵为 ⋀ m i {\bigwedge}_{m_i} ⋀mi​​的高斯分布，则单应性矩阵 H H H的极大似然估计可以通过最小化下面这个公式得到：

∑ i ( m i − m i ^ ) T ⋀ m i − 1 ( m i − m i ^ ) \sum_{i} (m_i-\hat{m_i})^T {\bigwedge}_{m_i}^{-1}(m_i-\hat{m_i}) i∑​(mi​−mi​^​)T⋀mi​−1​(mi​−mi​^​)

具体推导过程这里不再赘述。我们定义得到这个单应性矩阵为 H = [ h 1 h 2 h 3 ] H=\begin{bmatrix}h_1 & h_2 & h_3\end{bmatrix} H=[h1​​h2​​h3​​]，结合公式 ( 2 ) (2) (2)，可以得到：

[ h 1 h 2 h 3 ] = λ A [ r 1 r 2 t ] \begin{bmatrix} h_1 & h_2 & h_3 \end{bmatrix}= \lambda A \begin{bmatrix} r_1 & r_2 & t \end{bmatrix} [h1​​h2​​h3​​]=λA[r1​​r2​​t​]

即：

r 1 = 1 λ A − 1 h 1 r 2 = 1 λ A − 1 h 2 r_1=\frac{1}{\lambda}A^{-1}h_1 \\ r_2=\frac{1}{\lambda}A^{-1}h_2 r1​=λ1​A−1h1​r2​=λ1​A−1h2​

其中 λ \lambda λ是一个任意的标量。由于 r 1 r_1 r1​和 r 2 r_2 r2​是标准正交的，所以有：

r 1 T r 2 = h 1 T A − T A − 1 h 2 = 0 (3) r_1^T r_2=h_1^T A^{-T} A^{-1} h_2=0 \tag{3} r1T​r2​=h1T​A−TA−1h2​=0(3)

r 1 T r 1 = r 2 T r 2 = 1 ⇒ h 1 T A − T A − 1 h 1 = h 2 T A − T A − 1 h 2 (4) \begin{array}{cc} r_1^T r_1=r_2^T r_2=1 \\ \Rightarrow h_1^T A^{-T} A^{-1} h_1=h_2^{T} A^{-T} A^{-1} h_2 \end{array} \tag{4} r1T​r1​=r2T​r2​=1⇒h1T​A−TA−1h1​=h2T​A−TA−1h2​​(4)

给定一个单应性矩阵，对内参有上述两个基本约束条件。 由于单应性矩阵具有8个自由度和6个外参（3个旋转参数，3个平移参数），因此，到这我们只能得到内参的2个约束条件。

2 相机标定

本节详细介绍了如何有效解决相机标定问题。从一个解析解开始，然后是一个基于极大似然估计（可参考我的另一篇博文：极大似然估计）的非线性优化技术。最后，我们考虑了透镜的畸变，给出了解析解和非线性解。

2.1 封闭解

令：

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 − c α 2 β c v 0 − u 0 β α 2 β − c α 2 β c 2 α 2 β 2 + 1 β 2 − c ( c v 0 − u 0 β ) α 2 β 2 − v 0 β 2 c v 0 − u 0 β α 2 β − c ( c v 0 − u 0 β ) α 2 β 2 − v 0 β 2 ( c v 0 − u 0 β ) 2 α 2 β 2 + v 0 2 β 2 + 1 ] (5) \begin{aligned} 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{c}{\alpha^2\beta} & \frac{cv_0-u_0\beta}{\alpha^2\beta} \\ -\frac{c}{\alpha^2\beta} & \frac{c^2}{\alpha^2\beta^2}+\frac{1}{\beta^2} & -\frac{c(cv_0-u_0\beta)}{\alpha^2\beta^2}-\frac{v_0}{\beta^2} \\ \frac{cv_0-u_0\beta}{\alpha^2\beta} & -\frac{c(cv_0-u_0\beta)}{\alpha^2\beta^2}-\frac{v_0}{\beta^2} & \frac{(cv_0-u_0\beta)^2}{\alpha^2\beta^2}+\frac{v_0^2}{\beta^2}+1 \end{bmatrix} \end{aligned} \tag{5} B​=A−TA−1≡⎣⎡​B11​B12​B13​​B12​B22​B23​​B13​B23​B33​​⎦⎤​=⎣⎢⎡​α21​−α2βc​α2βcv0​−u0​β​​−α2βc​α2β2c2​+β21​−α2β2c(cv0​−u0​β)​−β2v0​​​α2βcv0​−u0​β​−α2β2c(cv0​−u0​β)​−β2v0​​α2β2(cv0​−u0​β)2​+β2v02​​+1​⎦⎥⎤​​(5)

B B B为对称矩阵，由一个6维的向量定义：

b = [ B 11 , B 12 , B 22 , B 13 , B 23 , B 33 ] T (6) b= \begin{bmatrix} B_{11},B_{12},B_{22},B_{13},B_{23},B_{33} \end{bmatrix} ^T \tag{6} b=[B11​,B12​,B22​,B13​,B23​,B33​​]T(6)

令矩阵 H H H中的第 i i i列向量为 h i = [ h i 1 , h i 2 , h i 3 ] T h_i=\begin{bmatrix}h_{i1},h_{i2},h_{i3}\end{bmatrix}^T hi​=[hi1​,hi2​,hi3​​]T，有：

h i T B h j = v i j T b (7) h_i^T B h_j=v_{ij}^T b \tag{7} hiT​Bhj​=vijT​b(7)

其中， v i j v_{ij} vij​为：

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 ] T v_{ij}= \begin{bmatrix} 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} \end{bmatrix} ^T vij​=[hi1​hj1​,hi1​hj2​+hi2​hj1​,hi2​hj2​,hi3​hj1​+hi1​hj3​,hi3​hj2​+hi2​hj3​,hi3​hj3​​]T

所以公式 ( 3 ) (3) (3)和 ( 4 ) (4) (4)就可以写为齐次形式：

[ v 12 T ( v 11 − v 22 ) T ] b = 0 (8) \begin{bmatrix} v_{12}^T \\ (v_{11}-v_{22})^T \end{bmatrix} b=0 \tag{8} [v12T​(v11​−v22​)T​]b=0(8)

如果获取 n n n张图像，也就相当于有 n n n个这样的方程，通过叠加这写方程，有：

V b = 0 (9) Vb=0 \tag{9} Vb=0(9)

其中V是一个2n×6的矩阵。如果n≥3，我们通常会有一个唯一解b。如果n=2，可以令倾斜参数c=0，即[0, 1, 0, 0, 0, 0]b = 0，这可作为额外的方程代入式(9)。

一旦估计出b，我们就可以计算出相机的内参矩阵A。有了内参矩阵，外参矩阵很容易计算。

2.2 极大似然估计

上面的解是通过最小化一个物理意义不大的代数距离得到的。我们可以通过极大似然估计来对它进行优化。

假设我们有棋盘格模型的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} \|m_{ij}-\hat{m}(A,R_i,t_i,M_j)\|^2 \tag{10} i=1∑n​j=1∑m​∥mij​−m^(A,Ri​,ti​,Mj​)∥2(10)

这里可以使用Levenberg-Marquardt算法进行迭代来求最优解，这个算法需要一个初始值，初始值可以从前面的描述中获取。

2.3 径向畸变校正

令 ( u , v ) (u,v) (u,v)为理想的无畸变的像素坐标， ( u ~ , v ~ ) (\tilde{u},\tilde{v}) (u~,v~)为对应的实际畸变后的像素坐标。 ( x , y ) (x,y) (x,y)为理想的无畸变的图像坐标， x ~ , y ~ \tilde{x},\tilde{y} x~,y~​为对应的畸变后的图像坐标。

实际情况中，可以用主点周围的泰勒级数展开的前几项进行描述。作者只考虑了径向畸变的前两项，因此有：

x ~ = x + x [ k 1 ( x 2 + y 2 ) + k 2 ( x 2 + y 2 ) 2 ] y ~ = y + y [ k 1 ( x 2 + y 2 ) + k 2 ( x 2 + y 2 ) 2 ] \begin{array}{cl} \tilde{x}=x+x[k_1(x^2+y^2)+k_2(x^2+y^2)^2] \\ \tilde{y}=y+y[k_1(x^2+y^2)+k_2(x^2+y^2)^2] \end{array} x~=x+x[k1​(x2+y2)+k2​(x2+y2)2]y~​=y+y[k1​(x2+y2)+k2​(x2+y2)2]​

其中 k 1 k_1 k1​和 k 2 k_2 k2​为径向畸变系数。径向畸变的中心与主点 ( u 0 , v 0 ) (u_0,v_0) (u0​,v0​)相同，因此有（主点周围的泰勒级数）：

u ~ = u + ( u − u 0 ) [ k 1 ( x 2 + y 2 ) + k 2 ( x 2 + y 2 ) 2 ] (11) \tilde{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) \tilde{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)

转为矩阵形式：

[ ( 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} \tilde{u}-u \\ \tilde{v}-v \end{bmatrix} [(u−u0​)(x2+y2)(v−v0​)(x2+y2)​(u−u0​)(x2+y2)2(v−v0​)(x2+y2)2​][k1​k2​​]=[u~−uv~−v​]

给定每张含有m个点的n张图像，我们可以联立所有的 2 m n 2mn 2mn个方程或写为矩阵形式： D k = d Dk=d Dk=d，其中 k = [ k 1 , k 2 ] T k=[k_1,k_2]^T k=[k1​,k2​]T。则线性最小二乘解可由下式给出：

k = ( D T D ) − 1 D T d (13) k=(D^TD)^{-1}D^Td \tag{13} k=(DTD)−1DTd(13)

扩展公式 ( 10 ) (10) (10)依然使用极大似然估计得到结果，利用LM优化算法计算使下列函数最小的参数值：

∑ i = 1 n ∑ j = 1 m ∥ m i j − m ^ ( A , k 1 , k 2 , R i , t i , M j ) ∥ 2 (14) \sum_{i=1}^{n} \sum_{j=1}^{m} \|m_{ij}-\hat{m}(A,k_1,k_2,R_i,t_i,M_j)\|^2 \tag{14} i=1∑n​j=1∑m​∥mij​−m^(A,k1​,k2​,Ri​,ti​,Mj​)∥2(14)