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

àì夳堔傛蜴生んèń 2023-07-06 13:04 25阅读 0赞

#### >>>文章索引<<< ####

[相机标定（一）-原理及内参、外参][-]  
相机标定（二）-畸变校正，张正友标定法  
[相机标定（三）-相机成像模型][- 1]

--------------------

# 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\[Rt\]M~, 其中A=⎣⎡α00cβ0u0v01⎦⎤(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描述了图像两个坐标轴之间的偏度。（关于相机内、外参矩阵的推导可参考：[链接][Link 1]）

本文中使用 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\[r1r2r3t\]⎣⎢⎢⎡XY01⎦⎥⎥⎤=A\[r1r2t\]⎣⎡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\[r1r2t\](2)

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

## 1.3 内参的约束条件 ##

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

假设 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=\[h1h2h3\]，结合公式 ( 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\} \[h1h2h3\]=λA\[r1r2t\]

即：

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=λ1A−1h1r2=λ1A−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\} r1Tr2=h1TA−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\} r1Tr1=r2Tr2=1⇒h1TA−TA−1h1=h2TA−TA−1h2(4)

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

# 2 相机标定 #

本节详细介绍了如何有效解决相机标定问题。从一个解析解开始，然后是一个基于极大似然估计（可参考我的另一篇博文：[极大似然估计][Link 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≡⎣⎡B11B12B13B12B22B23B13B23B33⎦⎤=⎣⎢⎡α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\} hiTBhj=vijTb(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=\[hi1hj1,hi1hj2\+hi2hj1,hi2hj2,hi3hj1\+hi1hj3,hi3hj2\+hi2hj3,hi3hj3\]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∑nj=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\]\[k1k2\]=\[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∑nj=1∑m∥mij−m^(A,k1,k2,Ri,ti,Mj)∥2(14)

[-]: https://blog.csdn.net/cfan927/article/details/104333527
[- 1]: https://blog.csdn.net/cfan927/article/details/105269221
[Link 1]: https://editor.csdn.net/md/?articleId=104333527
[Link 2]: https://blog.csdn.net/cfan927/article/details/104414285