3.相机标定原理（不考虑径向畸变）

3.相机标定原理（不考虑径向畸变）

3.1坐标转换

在进行相机标定之前，首先要进行相机的各个坐标的转换，最终得到得到投影矩阵。

1.相机标定:坐标系转换学习中已经详细的讲述了一个点从世界坐标系到像素坐标系下的坐标转换。

已知一点在世界坐标系中的坐标为$(x_w,y_w,z_w)$,其转换到像素坐标系下坐标$(P_{x_i},P_{x_i})$的转换公式如下：
$$\left[\begin{array}{c}{P}_{x_i}\ast Z_c\\ P_{y_i}\ast Z_c\\ Z_c\end{array}\right]=\left[\begin{array}{ccc}m\ast f& -m\ast f\ast cot\theta & m\ast P_x\\ 0& n\ast f/sin\theta & n\ast P_y\\ 0& 0& 1\end{array}\right]\ast (R|-Rc)\ast \left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]$$
这里所得到的坐标是$(P_{x_i},P_{x_i})$​的齐次坐标。

内参矩阵为了方便表述做如下修改:
$$K=\left[\begin{array}{ccc}m\ast f& -m\ast f\ast cot\theta & m\ast P_x\\ 0& n\ast f/sin\theta & n\ast P_y\\ 0& 0& 1\end{array}\right]====>\left[\begin{array}{ccc}\alpha & -\alpha \ast cot\theta & u_0\\ 0& \beta /sin\theta & v_0\\ 0& 0& 1\end{array}\right]$$
外参矩阵,其中$r_1,r_2,r_3$​为三维列向量：
$$(R|-Rc)====>[R|T]====>\left[\begin{array}{cc}{r}_1^T& t_x\\ r_2^T& t_y\\ r_3^T& t_z\end{array}\right]$$
投影矩阵M,其中$m_1,m_2,m_3$为四维行向量：
$$M=K\ast [R|t]=\left[\begin{array}{ccc}m\ast f& -m\ast f\ast cot\theta & m\ast P_x\\ 0& n\ast f/sin\theta & n\ast P_y\\ 0& 0& 1\end{array}\right]\ast (R|-Rc)=====>\left[\begin{array}{c}{m}_1\\ m_2\\ m_3\end{array}\right]$$

3.2标定目标

相机标定需要标定什么呢？

我们在得到上述的坐标转换公式之后，所需要做的就是确定出合理的相机内外参数，这样我们就可以利用该公式，通过世界坐标得到像素坐标了。

3.3标定过程

3.3.1标定装置

我们需要一个标定装置，（类似于棋盘格的东西），大概是这个样子的：

这个标定装置，我们以$O_w$​为世界坐标系的原点，每个格子的大小均已知，故在该坐标系下的每一点的位置坐标均已知。

3.3.2拍照

我们使用相机取拍摄这个标定装置，得到如上边右图的图片。

我们可以在图像中找到世界坐标系中一点的对应点，这个点在图像中的位置，即在像素坐标系中的位置是已知的。

3.3.3列方程组

为了方便区分，我们把世界坐标系中的点称为$P$​​，而把像素坐标系中的点称为$I$​​ 。由于我们根据坐标转换公式得到的坐标是$I$​​​​​的齐次坐标，这里我们要将其转换成欧式坐标，三维齐次坐标转化为二维欧式坐标，只需要除以齐次坐标的第三维就好了，举个例子：
$$\left[\begin{array}{c}{P}_{x_i}\ast Z_c\\ P_{y_i}\ast Z_c\\ Z_c\end{array}\right]=\left[\begin{array}{c}{m}_1\\ m_2\\ m_3\end{array}\right]\ast \left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right]=\left[\begin{array}{c}{m}_1\\ m_2\\ m_3\end{array}\right]\ast W=\left[\begin{array}{c}{m}_1\ast w\\ m_2\ast w\\ m_3\ast w\end{array}\right]$$
要将该坐标转换成欧式坐标：
$$\left[\begin{array}{c}{P}_{x_i}\\ P_{y_i}\end{array}\right]=\left[\begin{array}{c}(m_1\ast w)/(m_3\ast w)\\ (m_2\ast w)/(m_3\ast w)\end{array}\right]$$
同理，我们可以得到$I$​的欧式坐标,这里P为4×1的列向量(是世界坐标的齐次坐标)：
$$I=\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}(m_1\ast P)/(m_3\ast P)\\ (m_2\ast P)/(m_3\ast P)\end{array}\right]$$
于是我们就得到了两个方程：
$$\{\begin{array}{c}u=(m_1\ast P)/(m_3\ast P)---->-u\ast (m_3\ast P)+(m_1\ast P)=0\\ v=(m_2\ast P)/(m_3\ast P)---->-v\ast (m_3\ast P)+(m_2\ast P)=0\end{array}$$
我们已知M矩阵是一个3×4的矩阵，包含12个未知量，这十二个未知量由内参矩阵和外参矩阵中的参数共同组成。
内参矩阵包含有5个未知参数，外参矩阵包含12个位置参数，一共17个未知参数。但是外参矩阵有包含了6个约束条件，即旋转矩阵的3个列向量的模为1，且列向量互相正交。于是仅需要11个方程便可以求解得到M，即M矩阵的自由度为11。
因此我们至少需要6个点，才可以求解参数。而我们在实际的应用中通常会选取6个以上的点，来获得更加鲁棒(抗干扰)的结果。同时呢，我们还要保证这些点不能够全部来自同一个平面，否则会出现退化解。(不用纠结什么是退化解，注意不要让所有点来自同一个平面即可)
$$\begin{gathered} \{\begin{array}{c}-u_1\ast (m_3\ast P_1)+(m_1\ast P_1)=0\\ -v_1\ast (m_3\ast P_1)+(m_2\ast P_1)=0\end{array} \\ \{\begin{array}{c}-u_2\ast (m_3\ast P_2)+(m_1\ast P_2)=0\\ -v_2\ast (m_3\ast P_2)+(m_2\ast P_2)=0\end{array} \\ ... \\ \{\begin{array}{c}-u_n\ast (m_3\ast P_n)+(m_1\ast P_n)=0\\ -v_n\ast (m_3\ast P_n)+(m_2\ast P_n)=0\end{array} \end{gathered}$$
于是乎，我们就将相机标定的问题转化成了方程组求解的问题，只需要解出这个方程，就完成了标定。

3.3.4求解方程组，得到投影矩阵$M$

我们考虑将以上的方程组转换成矩阵方程。

这里我们需要定义两个矩阵$P,m$，这里的$P_n$为4×1的列向量,$m_1,m_2,m_3$为四维行向量
$$P={\left[\begin{array}{ccc}{P}_1^T& 0^T& -u_1{P}_1^T\\ 0& P_1^T& -V_1{P}_1^T\\ ...& ...& ...\\ P_n^T& 0^T& -u_n{P}_n^T\\ 0& P_n^T& -v_1{P}_n^T\end{array}\right]}_{2n\ast 12}m={\left[\begin{array}{c}{m}_1^T\\ m_2^T\\ m_3^T\end{array}\right]}_{12\ast 1}$$
于是我们就把方程组转换成了矩阵方程：$Pm=0$​​

已知上述方程组中的未知参数为11个，但方程的数量大于11个，于是就出现了超定方程组的问题，这导致我们的方程组无解。

于是我们的目标就转换为了寻找一个$m$使得$Pm$的值达到最小，即$min||Pm||$​。

我们容易发现，当$m$进行等比例缩小时，$Pm$的值会随着$m$的值的缩小而缩小，理论上当$m$趋向于无穷小时，$Pm$的值也会趋向于无穷小。显然，这不是我们想要的，我们需要一个确切的$m$的值。

因此，在这里我们还要引入一个约束条件$||m||=1$，这样我们就可以求取到确切的$m$​​的值了。

具体的求解方法是对$P$​进行奇异值分解，（奇异值分解？？矩阵论杀我！！这里对奇异值分解不做赘述，推导放在下边）

得出结论：$m$为$P$矩阵最小奇异值的右奇异向量，且$||m||=1$​

我们在求出$m$后，就可以得到$m_1,m_2,m_3$，于是乎我们就求解出了投影矩阵$M$​。这里将投影矩阵表示如下$a_1,a_2,a_3$为三维列向量：
$$M=\left[\begin{array}{c}{m}_1\\ m_2\\ m_3\end{array}\right]=[A|b]=\left[\begin{array}{cc}{a}_1^T& b_1\\ a_2^T& b_2\\ a_3^T& b_3\end{array}\right]=K\ast [R|T]$$

3.3.5求解内参数，外参数

上文中我们已经得到了投影矩阵$M$,根据$M=K[R|T]$可以求解内参数和外参数。推导过程复杂，先上结果吧
$$K=\left[\begin{array}{ccc}\alpha & -\alpha \ast cot\theta & u_0\\ 0& \beta /sin\theta & v_0\\ 0& 0& 1\end{array}\right],[R|T]=\left[\begin{array}{cc}{r}_1^T& t_x\\ r_2^T& t_y\\ r_3^T& t_z\end{array}\right]$$
由于我们所求到的$M$并非理想情况下的投影矩阵，因此在转换时需要一个转换系数ρ，来确保转换精度。
$$\begin{gathered} \rho M=K\ast [R|T] \\ \rho =\pm 1/|a_3| \end{gathered}$$
内参数：
$$\begin{gathered} \alpha ={\rho }^2|a_1\times a_3|sin\theta \\ \beta ={\rho }^2|a_2\times a_3|sin\theta \\ cos\theta =-(a_1\times a_3)\cdot (a_2\times a_3)/[|a_1\times a_3|\cdot |a_2\times a_3|] \\ {u}_0={\rho }^2(a_1\cdot a_3) \\ {v}_0={\rho }^2(a_2\cdot a_3) \end{gathered}$$
外参数：
$$\begin{gathered} r_1=(a_2\times a_3)/|a_2\times a_3| \\ {r}_3=\pm a_3/|a_3| \\ {r}_2=r_3\times r_1 \\ T=\rho K^{-1}b \end{gathered}$$

3.3.6推导过程

如下是数学公式的推导部分。(写到这里我已经后悔了o(╥﹏╥)o，公式太多了，就大概说一下吧)
在3.3.4中有提到对$P$​​​​​的奇异值分解。

$$\begin{gathered} P_{2n\ast 12}=U_{2n\ast 2n}{D}_{2n\ast 12}{V}_{12\ast 12}^H \\ D=\left[\begin{array}{cc}△& 0\\ 0& 0\end{array}\right] \end{gathered}$$
式中：$U，V$为酉矩阵，$D$​为奇异值矩阵,$△$​为奇异值对应的对角矩阵​

奇异值的算法：$P^{H}P$​的特征值为${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3\ge ...\ge {\lambda }_K-1\ge {\lambda }_k=...={\lambda }_n=0$​，我们称${\sigma }_i=\sqrt{{\lambda }_i}$​为矩阵$P$​​的奇异值​。于是我们就得到了矩阵$D$​。
$$\begin{gathered} V^H{P}^{H}PV=\left[\begin{array}{cc}{△}^2& 0\\ 0& 0\end{array}\right],V=(V_1,V_2) \\ (V_1,V_2{)}^H{P}^{H}P(V_1,V_2)=\left[\begin{array}{cc}{V}_1^H{P}^{H}P{V}_1& V_1^H{P}^{H}P{V}_2\\ V_2^H{P}^{H}P{V}_1& V_2^H{P}^{H}P{V}_2\end{array}\right]=\left[\begin{array}{cc}{△}^2& 0\\ 0& 0\end{array}\right] \end{gathered}$$
然后求$(V_1,V_2)$​：
$$\begin{gathered} \{\begin{array}{c}{V}_1^H{P}^{H}P{V}_1=(P{V}_1{)}^{H}P{V}_1={△}^2\\ V_2^H{P}^{H}P{V}_2=(P{V}_2{)}^{H}P{V}_2=0\end{array} \\ P{V}_2=0 \\ {△}^{-1}{V}_1^H{P}^{H}P{V}_1{△}^{-1}=I \end{gathered}$$
之后求$(U_1,U_2)$：
$$\begin{gathered} U=(U_1,U_2) \\ {U}_1=P{V}_1{△}^{-1} \\ {U}_2^H{U}_1=U_2^{H}P{V}_1{△}^{-1}=0(酉矩阵性质) \end{gathered}$$
$m$为$P$矩阵最小奇异值的右奇异向量，且$||m||=1$,即$V$中最小奇异值对应的奇异向量。

3.3.5中求内外参数的推导
$$\rho \left[\begin{array}{cc}{a}_1^T& b_1\\ a_2^T& b_2\\ a_3^T& b_3\end{array}\right]=\rho M=K\ast [R|T]=\left[\begin{array}{cc}\alpha r_1^T-\alpha \ast cot\theta r_2^T+u_0{r}_3^T& \alpha t_x-\alpha \ast cot\theta \ast t_y+u_0{t}_z\\ \beta r_2^T/sin\theta +v_0{r}_3^T& \beta t_y/sin\theta +v_0{t}_z\\ r_3^T& t_z\end{array}\right]$$
我们可以列出如下方程组：
$$\{\begin{array}{c}\rho a_1^T=\alpha r_1^T-\alpha \ast cot\theta r_2^T+u_0{r}_3^T--->\rho a_1=\alpha r_1-\alpha \ast cot\theta r_2+u_0{r}_3\\ \rho a_2^T=\beta r_2^T/sin\theta +v_0{r}_3^T--->\rho a_2=\beta r_2/sin\theta +v_0{r}_3\\ \rho a_3^T=r_3^T--->\rho a_3=r_3\end{array}$$
我们做如下计算：
( ρ a 1 ) × ( ρ a 3 ) = ρ 2 ( a 1 × a 3 ) = α r 2 − α ∗ c o t θ r 1 − − − > α = ρ 2 ( a 1 × a 3 ) / { r 2 − c o t θ r 1 } ( ρ a 2 ) × ( ρ a 3 ) = ρ 2 ( a 2 × a 3 ) = β r 1 / s i n θ − − − > β = ρ 2 ( a 2 × a 3 ) s i n θ / r 1 (ρa_1)×(ρa_3)=ρ^2(a_1×a_3)=αr_2-α*cotθr_1--->α=ρ^2(a_1×a_3)/\{r_2-cotθr_1\}\\ (ρa_2)×(ρa_3)=ρ^2(a_2×a_3)=βr_1/sinθ--->β=ρ^2(a_2×a_3)sinθ/r_1 (ρa1​)×(ρa3​)=ρ2(a1​×a3​)=αr2​−α∗cotθr1​−−−>α=ρ2(a1​×a3​)/{r2​−cotθr1​}(ρa2​)×(ρa3​)=ρ2(a2​×a3​)=βr1​/sinθ−−−>β=ρ2(a2​×a3​)sinθ/r1​

$$\begin{gathered} |\rho a_1|\times |\rho a_3|={\rho }^2|a_1\times a_3|=|\alpha |/sin\theta \\ |\rho a_2|\times |\rho a_3|={\rho }^2|a_2\times a_3|=|\beta |/sin\theta \end{gathered}$$

通过上述式子我们可以的到：
$$\begin{gathered} \alpha ={\rho }^2|a_1\times a_3|sin\theta \\ \beta ={\rho }^2|a_2\times a_3|sin\theta \end{gathered}$$
做如下计算我们可以求得$u_0,v_0,cos\theta$都是通过点乘和叉乘得到的（就不推导了吧）：
$$\begin{gathered} u_0={\rho }^2(a_1\cdot a_3) \\ {v}_0={\rho }^2(a_2\cdot a_3) \\ cos\theta =-(a_1\times a_3)\cdot (a_2\times a_3)/[|a_1\times a_3|\cdot |a_2\times a_3|] \end{gathered}$$
已知我们的旋转矩阵$R$​​​是一个正交矩阵，因此其行向量和列向量的模均为1，于是可以求解得到$\rho$：
$$\begin{gathered} \rho a_3=r_3 \\ \rho |a_3|=|r_3|=\pm 1 \\ \rho =\pm 1/|a_3| \end{gathered}$$
关于外参数，我们在上述求$\alpha ,\beta$​​​的推导中已经求得了$r_1,r_3$​：
$$\begin{gathered} r_1=(a_2\times a_3)/|a_2\times a_3| \\ {r}_3=\pm a_3/|a_3| \end{gathered}$$
同样利用旋转矩阵是正交矩阵的性质，可以得到：
$$r_2=r_3\times r_1$$
由 $\rho b=KT$​ 易知，$T=K^{-1}\rho b=\rho K^{-1}b$​​，由于$K,b$均为已知的，可求得$T$​

推导完毕！

$$\begin{gathered} u_0={\rho }^2(a_1\cdot a_3) \\ {v}_0={\rho }^2(a_2\cdot a_3) \\ cos\theta =-(a_1\times a_3)\cdot (a_2\times a_3)/[|a_1\times a_3|\cdot |a_2\times a_3|] \end{gathered}$$
已知我们的旋转矩阵$R$​​​是一个正交矩阵，因此其行向量和列向量的模均为1，于是可以求解得到$\rho$：
$$\begin{gathered} \rho a_3=r_3 \\ \rho |a_3|=|r_3|=\pm 1 \\ \rho =\pm 1/|a_3| \end{gathered}$$
关于外参数，我们在上述求$\alpha ,\beta$​​​的推导中已经求得了$r_1,r_3$​：
$$\begin{gathered} r_1=(a_2\times a_3)/|a_2\times a_3| \\ {r}_3=\pm a_3/|a_3| \end{gathered}$$
同样利用旋转矩阵是正交矩阵的性质，可以得到：
$$r_2=r_3\times r_1$$
由 $\rho b=KT$​ 易知，$T=K^{-1}\rho b=\rho K^{-1}b$​​，由于$K,b$均为已知的，可求得$T$​

推导完毕！

参考资料

计算机视觉 鲁鹏