Geometric Camera Calibration

Geometric camera calibration using circular control points

链接http://www.ee.oulu.fi/~jth/calibr/

作者：竹石

相机标定可以归纳为$P-n-P$​​​​(Perspective-n-Point)的问题，即已知三维物点坐标和对应的二维投影坐标，求解相机参数。由于镜头的畸变(径向和切向)带来非线性成像模型，一般求解方法分为两步：

• 不考虑畸变，成像模型为线性模型，利用线性求解方法求出初始解
• 考虑畸变，利用初始解和成像模型对三维物点投影得到的投影点与观测点形成最优问题， 通过最小二乘进行估计。

这篇文章的精彩之处在于给出逆畸变模型，在上两步的基础上，利用逆畸变模型进一步优化畸变参数。

文章的主要框架内容：

1.相机模型

1.1正投影模型

$$\left[\begin{array}{l}u\\ v\\ 1\end{array}\right]=\lambda P\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{l}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]=\lambda F\left[\begin{array}{l}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right]\text{\hspace{1em}(1)}$$

其中，$A_w=[X_w{Y}_w{Z}_w{]}^T$​​ 表示世界坐标系某一物点，$a=[u \quad v]^T $ 表示图像坐标系中的投影点，$\lambda$​表示缩放因子。

相机的内参：
$$P=\left[\begin{array}{llll}sf& 0& u_0& 0\\ 0& f& v_0& 0\\ 0& 0& 1& 0\end{array}\right]$$
其中$f$表示焦距，$(u_0,v_0)$表示主点坐标，$s$​表示横纵比。

相机的外参：

旋转矩阵 $R$​ 是一个$3\times 3$​的正交旋转矩阵，可用三个欧拉角$w$​、$\phi$​ 、$k$​ 定义；平移矩阵$t=[t_x{t}_y{t}_z]$​​。

相机的畸变模型：
$$a_c=a_d+{\mathcal{F}}_D\left({\mathbf{a}}_d,\delta \right)\text{\hspace{1em}(2)}$$
其中$a_c=[u_c{v}_c{]}^T$表示矫正后的投影点，$a_d=[u_d{v}_d{]}^T$表示畸变点。${\mathcal{F}}_D\left({\mathbf{a}}_d,\delta \right)$​为多项式，表示如下，
$${\mathcal{F}}_D\left({\mathbf{a}}_d,\delta \right)=\left[\begin{array}{l}{\bar{u}}_d\left(k_1{r}_d^2+k_2{r}_d^4+k_3{r}_d^6+\dots \right)+\left(2p_1{\bar{u}}_d{\bar{v}}_d+p_2\left(r_d^2+2{\bar{u}}_d^2\right)\right)\left(1+p_3{r}_d^2+\dots \right)\\ {\bar{v}}_d\left(k_1{r}_d^2+k_2{r}_d^4+k_3{r}_d^6+\dots \right)+\left(p_1\left(r_d^2+2{\bar{v}}_d^2\right)+2p_2{\bar{u}}_d{\bar{v}}_d\right)\left(1+p_3{r}_d^2+\dots \right)\end{array}\right]$$
其中$\delta ={\left[k_1,k_2,\dots ,p_1,p_2,\dots \right]}^{\mathrm{T}}$为畸变系数，$k_1，k_2,…… $表示径向畸变系数，$p_1，p_2,…… $表示切向畸变系数，$r_{c}^{2}=u_{c}{2}+v_{c}^{2}$。

1.2反投影模型

设平面$\Pi$​ 以$h_0=[X_0{Y}_0{Z}_0]$​为原点，以$h_1=[X_1{Y}_1{Z}_1]$​和$h_2=[X_2{Y}_2{Z}_2]$​为两个基。则平面$\Pi$​上的任意一点$(X_H,Y_H)$​​对应的世界坐标系为：$X_H{h}_1+Y_H{h}_2+h_0$​​。

即：
$$\left[\begin{array}{ccc}{h}_1& h_2& h_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{X}_H\\ Y_H\\ 1\end{array}\right]=H\left[\begin{array}{c}{X}_H\\ Y_H\\ 1\end{array}\right]$$
正投影模型为：
$$\left[\begin{array}{c}{u}_H\\ v_H\\ 1\end{array}\right]=F{P}_H=\lambda FH\left[\begin{array}{c}{X}_H\\ Y_H\\ 1\end{array}\right]$$
即反投影模型：
$$\lambda \left[\begin{array}{c}{X}_H\\ Y_H\\ 1\end{array}\right]=(FH{)}^{-1}\left[\begin{array}{c}{u}_H\\ v_H\\ 1\end{array}\right]\text{\hspace{1em}(3)}$$

1.3需要标定的参数：

6个相机外参：
$$t_x、t_y、t_z、w、\phi 、k$$
8个相机内参(仅考虑了径向和切向起主导的前两个畸变系数)：
$$s、f、u_0、v_0、k_1、k_2、p_1、p_2$$

2.圆形标定点的偏差校正

透视投影不是保形变换，直线在透视投影模型下为直线，一般二维或三维形状与图像平面不共面时会发生变形。常用的标定板是棋盘格，棋盘格的角点是包型变换，但不易精准检测。圆形标定板也是校准中常用的标定板，圆形可以准确的找到中心点，但通过透视投影圆心会发生偏差。

平面二次曲线一般形式:
$$A{X}_H^2+2B{X}_H{Y}_H+C{Y}_H^2+2D{X}_H+2E{Y}_H+F=0$$
矩阵形式:
$$\left[\begin{array}{lll}{X}_H& Y_H& 1\end{array}\right]\left[\begin{array}{ccc}A& B& D\\ B& C& E\\ D& E& F\end{array}\right]\left[\begin{array}{c}{X}_H\\ Y_H\\ 1\end{array}\right]=0$$
令：
$$Q=\left[\begin{array}{lll}A& B& D\\ B& C& E\\ D& E& F\end{array}\right]$$
则圆的表示形式：
$$\begin{array}{r}Q=\left[\begin{array}{ccc}-\frac{1}{r^2}& 0& 0\\ 0& -\frac{1}{r^2}& 0\\ 0& 0& 1\end{array}\right]\end{array}$$
得：
$$\left[\begin{array}{lll}{X}_H& Y_H& 1\end{array}\right]\mathrm{Q}\left[\begin{array}{c}{X}_H\\ Y_H\\ 1\end{array}\right]=0\text{\hspace{1em}(4)}$$
因为反相机模型：
$$\lambda \left[\begin{array}{c}{X}_H\\ Y_H\\ 1\end{array}\right]=(FH{)}^{-1}\left[\begin{array}{c}{u}_H\\ v_H\\ 1\end{array}\right]\text{\hspace{1em}(5)}$$
（5）带入（4）得

$${\left(\begin{array}{l}(FH{)}^{-1}\left[\begin{array}{c}{u}_H\\ v_H\\ 1\end{array}\right]\end{array}\right)}^{T}Q(FH{)}^{-1}\left[\begin{array}{c}{u}_H\\ v_H\\ 1\end{array}\right]=0$$
进一步整理:
$${\left(\left[\begin{array}{c}{u}_H\\ v_H\\ 1\end{array}\right]\right)}^T{\left((FH{)}^{-1}\right)}^{T}Q(FH{)}^{-1}\left[\begin{array}{c}{u}_H\\ v_H\\ 1\end{array}\right]=0$$
令$S={\left((FH{)}^{-1}\right)}^{T}Q(FH{)}^{-1}$，则$S^{-1}=FH{Q}^{-1}{H}^T{F}^T$​。​

令$G=H{Q}^{-1}{H}^T$，由引文[1]得：
$$\lambda \left[\begin{array}{l}{u}_e\\ v_e\\ 1\end{array}\right]=FG{F}^T\left[\begin{array}{l}0\\ 0\\ 1\end{array}\right]=FG{f}_3\text{\hspace{1em}(6)}$$
其中$f_3$ 表示 $F$ 的最后一行，利用(6)代替相机正投影模型(1)可以修正椭圆的中心，获得无偏估计。

3.逆畸变模型

畸变矫正：
$$a_c=a_d+{\mathcal{F}}_D\left(a_d,\delta \right)\text{\hspace{1em}(7)}$$
其中, 正确或矫正后的坐标 $a_c={\left[u_c,v_c\right]}^T;$ 畸变坐标: $a_d={\left[u_d,v_d\right]}^T$​ 。

3.1递归逆畸变模型

由(4)可得:
$$a_d=a_c-{\mathcal{F}}_D\left(a_d,\delta \right)=a_c-{\mathcal{F}}_D\left(a_c-{\mathcal{F}}_D\left(a_d,\delta \right),\delta \right)=\cdots$$
在每次迭代中，等号右侧用(7)替换$a_d$​​时会使误差减少。在实践中，至少需要三到四次迭代来补偿镜头畸变，需要在图像的不同位置多次评估失真函数，使得递归逆畸变模型不实用。

3.2非递归逆畸变模型:

采用一阶泰勒级数对${\mathcal{F}}_D\left(a_d,\delta \right)$​近似表示:
$$\begin{array}{rl}& {\mathcal{F}}_D\left(a_d,\delta \right)\approx {\mathcal{F}}_D\left(a_c,\delta \right)+D\left(a_c\right)\left(a_d-a_c\right)\end{array}$$
(7)可以表示为：
$$a_c\approx a_d+{\mathcal{F}}_D\left(a_c,\delta \right)+D\left(a_c\right)\left(a_d-a_c\right)$$
则，
$$a_d\approx a_c-\frac{1}{d_{11}\left(a_c\right)+d_{22}\left(a_c\right)+1}{\mathcal{F}}_D\left(a_c,\delta \right)\equiv a_c-{\mathcal{F}}_D^{\ast }\left(a_c,\delta \right)$$
两个径向 $k_1$ 和 $k_2$ 及切向畸变系数 $p_1$ 和 $p_2$, 则逆畸变模型的近似值表示:
$${\mathcal{F}}_D^{\ast }\left(a_c,\delta \right)\approx \frac{{\mathcal{F}}_D\left(a_c,\delta \right)}{4k_1{r}_c^2+6k_2{r}_c^4+8p_1{v}_c+8p_2{u}_c+1}\text{\hspace{1em}(8)}$$
其中, $r_c^2=u_c^2+v_c^2$
那么对像素点 $a_c$, 利用逆畸变模型进行畸变后的坐标表示为 $a_d$, 数学描述为：
$$a_d\approx a_c-{\mathcal{F}}_D^{\ast }\left(a_c,\delta \right)\text{\hspace{1em}(9)}$$

4.利用逆畸变模型优化畸变系数

通过step1和step2得到的相机参数，利用公式(7)进行矫正,得到矫正坐标$a_c$​​。利用(9)对矫正坐标$a_c$​进行扭曲得到扭曲坐标，则：

$$a_c-a_d={\mathcal{F}}_D\left(a_d,\delta \right)=B\delta \text{\hspace{1em}(10)}$$
畸变系数:
$${\mathcal{F}}_D\left({\mathbf{a}}_d,\delta \right)=\left[\begin{array}{l}{\bar{u}}_d\left(k_1{r}_d^2+k_2{r}_d^4+k_3{r}_d^6+\dots \right)+\left(2p_1{\bar{u}}_d{\bar{v}}_d+p_2\left(r_d^2+2{\bar{u}}_d^2\right)\right)\left(1+p_3{r}_d^2+\dots \right)\\ {\bar{v}}_d\left(k_1{r}_d^2+k_2{r}_d^4+k_3{r}_d^6+\dots \right)+\left(p_1\left(r_d^2+2{\bar{v}}_d^2\right)+2p_2{\bar{u}}_d{\bar{v}}_d\right)\left(1+p_3{r}_d^2+\dots \right)\end{array}\right]$$
$B$ 矩阵可以写为:
$$\mathbf{B}(i)=\left[\begin{array}{lllllc}{\bar{u}}_d(i)r_d^2(i)& {\bar{u}}_d(i)r_d^4(i)& \dots & 2{\bar{u}}_d(i){\bar{v}}_d(i)& r_d^2(i)+2{\bar{u}}_d^2(i)& \dots \\ {\bar{v}}_d(i)r_d^2(i)& {\bar{v}}_d(i)r_d^4(i)& \dots & r_d^2(i)+2{\bar{v}}_d^2(i)& 2{\bar{u}}_d(i){\bar{v}}_d(i)& \dots \end{array}\right]$$
$\delta$ 矩阵可以写为:
$$\delta =\left[\begin{array}{llllll}{k}_1& k_2& \dots & p_1& p_2& \dots \\ k_1& k_2& \dots & p_1& p_2& \dots \end{array}\right]$$
则畸变系数优化求解方法:
$$\delta =[B{]}^{+}\left[a_c-a_d\right]$$
其中，$[.{]}^{+}$表示矩阵的伪逆。

5.验证逆畸变模型的精度

求未矫正坐标和扭曲坐标的距离$|a_d-a_0|$​​​​来验证逆畸变模型的正确性。

未矫正坐标和扭曲坐标之间的差异可用直方图表示,表明误差小于0.01像素。

引文：

[1]K. Kanatani,Geometric computation for Machine Vision, Oxford: Clarendon Press, 1993.