【笔记】基于图像的三维重建——相机几何与相机标定

本文主要关于深蓝学院系列课程——基于图像的三维重建的笔记。

课程链接 基于图像的三维重建

1、相机几何

1.1、针孔模型 & 透镜

• 针孔模型：光线的传播导致 重叠，故需要遮挡过多的光线，使得 每一个点在胶片上只有一个成像点。

• 针孔相机的数学模型： 利用相似三角形的关系，可以得到真实空间中的一点与成像点的映射关系。

$$\boldsymbol{P}=\begin{pmatrix}x\\y\\z\end{pmatrix}\to\boldsymbol{p}=\begin{pmatrix}x'\\y'\end{pmatrix}\quad\begin{cases}x'=f\dfrac{x}{z}\\y'=f\dfrac{y}{z}\end{cases}$$

• 针孔大小对成像的影响： 光圈大，光线亮，成像模糊 （多条光线在同一成像点）；光圈小，光线暗，成像清晰。

• 透镜： 使得同一个点的多条光线可以聚焦到胶片上，增加照片亮度。

• 透镜问题：

1. 失焦：透镜自身工艺原因不能将光线很好地汇集。可以清晰成像的部分称为景深。

1. 径向畸变：光线在远离透镜中心的地方比靠近中心的地方更加弯曲。主要分为 桶型和枕型 畸变。

1.2、相机几何

• 相机几何模型涉及到的坐标系： 相机坐标系、成像坐标系、像素坐标系。

• 成像坐标系与像素坐标系的转换关系： 统一单位系数 $k, l$ ，坐标系偏置 $c_x, c_y$.

• 齐次坐标系的引入和投影矩阵： 上述变换是 非线性的，在表达上十分不方便（非线性变换），为使其可以 线性表示，引入齐次坐标系表示。$M$ 称为相机的投影矩阵，表示三维点投影到 像素坐标系 的变换关系。

$$p_h=\begin{pmatrix}\alpha x+c_xz\\\beta y+c_yz\\z\end{pmatrix}=\begin{pmatrix}\alpha&0&c_x&0\\0&\beta&c_y&0\\0&0&1&0\end{pmatrix}\boxed{\begin{pmatrix}x\\y\\z\\1\end{pmatrix}}—P_h=MP_h$$

• 相机偏斜： 像素坐标系不满足严格的正交关系。

$$\boldsymbol{p}=\begin{pmatrix}\alpha&-\alpha\cot\theta&c_x&0\\0&\beta/\sin\theta&c_y&0\\0&0&1&0\end{pmatrix}\begin{pmatrix}x\\y\\z\\1\end{pmatrix}$$

• 相机坐标系下的相机模型：

$$p=K(I\quad0)P=K(I\quad 0)\begin{pmatrix}R&t\\0^\mathrm{T}&1\end{pmatrix}P_\mathrm{w}=K(R\quad t)P_\mathrm{w}=MP_\mathrm{w}\\ 相机内参：\boldsymbol{K}=\begin{pmatrix}\alpha&&-\alpha\cot\theta&&c_x\\0&&\beta/\sin\theta&&c_y\\0&&0&&1\end{pmatrix}\\ 相机外参：({R}\quad t)$$

其中 $M$ 为投影矩阵； $K$ 为 相机内参，决定了相机坐标系下的空间点到图像点的映射； $({R}\quad t)$ 为 相机外参，表示相机坐标系相对世界坐标系的旋转和平移。

• 相机模型内外参的自由度： 相机模型一共 11个自由度，其中内参 5个自由度 ，外参 6个自由度（三维旋转+平移）。

• $p$ 转换为欧式坐标：

$${p=MP_w=\begin{pmatrix}m_1^\mathrm{T}\\m_2^\mathrm{T}\\m_3^\mathrm{T}\end{pmatrix}P_w=\begin{pmatrix}m_1^\mathrm{T}P_w\\m_2^\mathrm{T}P_w\\m_3^\mathrm{T}P_w\end{pmatrix}}\to\left(\frac{m_1^\mathrm{T}P_w}{m_3^\mathrm{T}P_w},\frac{m_2^\mathrm{T}P_w}{m_3^\mathrm{T}P_w}\right)^\mathrm{T}\\ M=\begin{pmatrix}m_1^\mathrm{T}\\m_2^\mathrm{T}\\m_3^\mathrm{T}\end{pmatrix}$$

• 相关定理：

$$M=K(R\quad t)=(KR\quad Kt)=(A\quad b)\quad A=\begin{pmatrix}a_1^\mathrm{T}\\a_2^\mathrm{T}\\a_3^\mathrm{T}\end{pmatrix}$$

1. $M$ 是透视投影矩阵的一个充分必要条件是 $Det(A)\ne 0$.
2. $M$ 是零倾斜透视投影矩阵的一个充分必要条件是 $Det(A)\ne 0$ 且 $(a_1\times a_3)\cdot(a_2\times a_3)=0$.
3. $M$ 是零倾斜且宽高比为 1 的 透视投影矩阵的一个充分必要条件是 $Det(A)\ne 0$ 且 $\left\{\begin{matrix}(a_1\times a_3)\cdot(a_2\times a_3)=0\(a_1\times a_3)\cdot(a_1\times a_3)=(a_2\times a_3)\cdot(a_2\times a_3)\end{matrix}\right.$

• 投影变换的性质：
  • 点投影为点；
  • 线投影为线；
  • 近大远小；
  • 角度不再保持，平行线相交。

1.3、其他相机模型

• 弱透视与透视投影相机： 弱透视是当物体离光心无限远，所有成像点 可以视为在同一平面。

• 思考： 缺少了哪三维的自由度？

  • 答：关于 $X，Y$ 轴的旋转以及 $Z$ 轴的平移。
• 思考： 最后一行为何是 $(0,0,0,1)$？
  • 参考：三维重建1——相机几何模型和投影矩阵

• 正交投影：光心距离成像平面以及物体均无限远。应用与工业设计，建筑设计。

2、相机标定

2.1、针孔模型 & 透镜相机标定问题

目标： 从1张或者多张图像中估算内、外参数矩阵。

方法： 通过标定装置，给定世界坐标系下的坐标点，同时从图像上得到坐标点在像素坐标系下的坐标，求解投影矩阵。

已知： 待求参数量（未知量）为 11个，一对对应点可以提供 2个方程。

求解： 至少需要 6对点对 进行求解。实际操作一般采用多于6对点对来获得更为鲁棒的结果。

$$\begin{aligned}u_i&=\frac{\boldsymbol{m}_1^\mathrm{T}\boldsymbol{P}_i}{\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i}\quad\to\quad u_i\big(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i\big)=\boldsymbol{m}_1^\mathrm{T}\boldsymbol{P}_i\quad\to\quad\boldsymbol{m}_1^\mathrm{T}\boldsymbol{P}_i-u_i(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i)=0\\\\v_i&=\frac{\boldsymbol{m}_2^\mathrm{T}\boldsymbol{P}_i}{\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i}\quad\to\quad v_i\big(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i\big)=\boldsymbol{m}_2^\mathrm{T}\boldsymbol{P}_i\quad\to\quad\boldsymbol{m}_2^\mathrm{T}\boldsymbol{P}_i-v_i(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i)=0\end{aligned}$$

标定问题： 方程数 $2n$ 个且 $n\ge 6$ ，未知参数 11个。为 超定齐次线性方程组 。

$$\boldsymbol{P}\stackrel{\mathrm{def}}{=}\begin{pmatrix}\boldsymbol{P}_{1}^{\mathrm{T}}&\boldsymbol{0}^{\mathrm{T}}&-u_{1}\boldsymbol{P}_{1}^{\mathrm{T}}\\\boldsymbol{0}^{\mathrm{T}}&\boldsymbol{P}_{1}^{\mathrm{T}}&-v_{1}\boldsymbol{P}_{1}^{\mathrm{T}}\\\cdots&\cdots&\cdots\\\boldsymbol{P}_{n}^{\mathrm{T}}&\boldsymbol{0}^{\mathrm{T}}&-u_{n}\boldsymbol{P}_{n}^{\mathrm{T}}\\\boldsymbol{0}^{\mathrm{T}}&\boldsymbol{P}_{n}^{\mathrm{T}}&-v_{n}\boldsymbol{P}_{n}^{\mathrm{T}}\end{pmatrix}_{2n\times12},\boldsymbol{m}\triangleq\begin{pmatrix}\boldsymbol{m_1}\\\boldsymbol{m_2}\\\boldsymbol{m_3}\end{pmatrix}_{12\times1} \\ \begin{cases}-u_1(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_1)+\boldsymbol{m}_1^\mathrm{T}\boldsymbol{P}_1=0\\-v_1(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_1)+\boldsymbol{m}_2^\mathrm{T}\boldsymbol{P}_1=0\\\vdots\\-u_n(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_n)+\boldsymbol{m}_1^\mathrm{T}\boldsymbol{P}_n=0\\-v_n(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_n)+\boldsymbol{m}_2^\mathrm{T}\boldsymbol{P}_n=0\end{cases} \to Pm=0$$

转换为最小二乘问题：

$$\begin{aligned}Pm=0\quad\Longleftrightarrow\quad&\min_m\|Pm\|\\&s.t.\|m\|=1\end{aligned}$$

求解上式：

1. 矩阵 $P$ 进行奇异值分解 $P=UDV^T$
2. $x^*$ 为$V$ 矩阵的最后一列（最小特征值对应的右特征向量）

求解参数： 式中 $\rho$ 为 $\left | m \right | = 1$ 的约束带来的。正是由于这个约束，使得 $M$ 矩阵为11个自由度即可！

$$\rho \boldsymbol{M} =\boldsymbol{K}(\boldsymbol{R}\quad\boldsymbol{t})=\begin{pmatrix}\alpha\boldsymbol{r}_{1}^{\mathrm{T}}-\alpha\cot\theta\boldsymbol{r}_{2}^{\mathrm{T}}+c_{x}\boldsymbol{r}_{3}^{\mathrm{T}}&\alpha t_{x}-\alpha\cot\theta t_{y}+c_{x}t_{z}\\\frac{\beta}{\sin\theta}\boldsymbol{r}_{2}^{\mathrm{T}}+c_{y}\boldsymbol{r}_{3}^{\mathrm{T}}&\frac{\beta}{\sin\theta}t_{y}+c_{y}t_{z}\\\boldsymbol{r}_{3}^{\mathrm{T}}&t_{z}\end{pmatrix}$$

前置变换：

$$\rho(A\quad b)=K(R\quad t)\\ A=\begin{pmatrix}\boldsymbol{a}_{1}^{\mathrm{T}}\\\boldsymbol{a}_{2}^{\mathrm{T}}\\\boldsymbol{a}_{3}^{\mathrm{T}}\end{pmatrix}\quad\boldsymbol{b}=\begin{pmatrix}b_{1}\\b_{2}\\b_{3}\end{pmatrix}$$

求得结果：

$$\rho=\frac{\pm1}{|a_{3}|}\\\ c_{x}=\rho^{2}(a_{1}\cdot a_{3}) \\ c_{y}=\rho^{2}(a_{2}\cdot a_{3}) \\ \cos\theta=-{\frac{(a_{1}\times a_{3})\cdot(a_{2}\times a_{3})}{|a_{1}\times a_{3}|\cdot|a_{2}\times a_{3}|}} \\ \alpha=\rho^{2}|\boldsymbol{a}_{1}\times\boldsymbol{a}_{3}|\sin\theta \\ \beta=\rho^2|\boldsymbol{a}_2\times\boldsymbol{a}_3|\sin\theta \\ r_{1}={\frac{(a_{2}\times a_{3})}{|a_{2}\times a_{3}|}} \\ r_{2}=r_{3}\times r_{1} \\ r_{3}={\frac{\pm a_{3}}{|a_{3}|}}\\ \boldsymbol{t}=\rho\boldsymbol{K}^{-1}\boldsymbol{b}$$

退化情况： 所选取点对不能位于同一平面。

2.2、径向畸变的相机标定

建模： 畸变示意图如下

建模方程如下： 其中 $k_p$ 为畸变因子

$$\\ \begin{pmatrix}\dfrac{1}{\lambda}&0&0\\0&\dfrac{1}{\lambda}&0\\0&0&1\end{pmatrix}M\boldsymbol{P}_i\to\binom{u_i}{v_i}=\boldsymbol{p}_i \\ \lambda=1\pm\sum_{p=1}^{q}k_{p}d^{2p}，\quad d^2=u^2+v^2$$

设畸变参数矩阵与 $M$ 为 $Q$ ：

$${\boxed{\begin{pmatrix}\frac1\lambda&0&0\\0&\frac1\lambda&0\\0&0&1\end{pmatrix}M}}P_i\to\binom{\lambda_i}{\nu_i}=p_i\quad Q=\begin{pmatrix}q_1^\mathrm{T}\\q_2^\mathrm{T}\\q_3^\mathrm{T}\end{pmatrix}\\$$

构建方程组，该方程组为 非线性 方程组：

$$p_i=\begin{pmatrix}u_i\\v_i\end{pmatrix}=\begin{pmatrix}\dfrac{q_1^\mathrm{T}P_i}{q_3^\mathrm{T}P_i}\\\dfrac{q_2^\mathrm{T}P_i}{q_3^\mathrm{T}P_i}\end{pmatrix}\longrightarrow\begin{cases}u_iq_3^\mathrm{T}P_i=q_1^\mathrm{T}P_i\\v_iq_3^\mathrm{T}P_i=q_2^\mathrm{T}P_i\end{cases}$$

方法一： 采用 牛顿法或L-M法 进行非线性求解：

$$\begin{cases}u_iq_3^\mathrm{T}P_i=q_1^\mathrm{T}P_i\\v_iq_3^\mathrm{T}P_i=q_2^TP_i\end{cases}\longrightarrow\begin{cases}u_iq_3^\mathrm{T}P_i-q_1^\mathrm{T}P_i=0\\v_iq_3^\mathrm{T}P_i-q_2^\mathrm{T}P_i=0\end{cases}\longrightarrow\begin{cases}f_1(k_1,k_2,k_3,m_1^\mathrm{T},m_2^\mathrm{T},m_3^\mathrm{T})=0\\f_2(k_1,k_2,k_3,m_1^\mathrm{T},m_2^\mathrm{T},m_3^\mathrm{T})=0\end{cases}\\$$ $$\begin{array}{l}\text{令}x^\mathrm{T}=(k_1,k_2,k_3,\boldsymbol{m}_1^\mathrm{T},\boldsymbol{m}_2^\mathrm{T},\boldsymbol{m}_3^\mathrm{T})\\\\x^*=arg\min_x\|f(x)\|^2\end{array}$$

方法二： 先求解线性部分，再求解非线性部分：

$$\boldsymbol{p}_i=\binom{u_i}{v_i}=\frac1\lambda\begin{pmatrix}\cfrac{\boldsymbol{m}_1^\mathrm{T}\boldsymbol{P}_i}{\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i}\\\cfrac{\boldsymbol{m}_2^\mathrm{T}\boldsymbol{P}_i}{\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i}\end{pmatrix}\longrightarrow\quad\frac{u_i}{v_i}=\frac{\frac1\lambda\frac{\left(\boldsymbol{m}_1^\mathrm{T}\boldsymbol{P}_i\right)}{\left(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i\right)}}{\frac1\lambda\frac{\left(\boldsymbol{m}_2^\mathrm{T}\boldsymbol{P}_i\right)}{\left(\boldsymbol{m}_3^\mathrm{T}\boldsymbol{P}_i\right)}}=\frac{\boldsymbol{m}_1^\mathrm{T}\boldsymbol{P}_i}{\boldsymbol{m}_2^\mathrm{T}\boldsymbol{P}_i}$$ $$\begin{cases}v_1(\boldsymbol{m_1^\mathrm{T}}\boldsymbol{P_1})-u_1(\boldsymbol{m_2^\mathrm{T}}\boldsymbol{P_1})=0\\\vdots\\v_i(\boldsymbol{m_1^\mathrm{T}}\boldsymbol{P_i})-u_i(\boldsymbol{m_2^\mathrm{T}}\boldsymbol{P_i})=0\\\vdots\\v_n(\boldsymbol{m_1^\mathrm{T}}\boldsymbol{P_n})-u_n(\boldsymbol{m_2^\mathrm{T}}\boldsymbol{P_n})=0\end{cases} ,\quad L\stackrel{\mathrm{def}}{=}\begin{pmatrix}v_1\boldsymbol{P_1^\mathrm{T}}&-u_1\boldsymbol{P_1^\mathrm{T}}\\v_2\boldsymbol{P_2^\mathrm{T}}&-u_2\boldsymbol{P_2^\mathrm{T}}\\\vdots&\vdots\\v_n\boldsymbol{P_n^\mathrm{T}}&-u_n\boldsymbol{P_n^\mathrm{T}}\end{pmatrix}$$ $$构建 \quad Ln = 0，通过SVD求得 \quad n=\binom{m_1}{m_2}$$

非线性部分：

$m_3,k_1,k_2,k_3$ 关于 $m_1,m_2$ 的非线性函数，采用 牛顿法和L-M法 求解。