重拾图形图像处理 ---- 笔试面试题：三维重建相关（1）

文章大纲

齐次坐标、点到直线距离
给三角形三边求面积
简述SIFT特征点检测、描述和匹配的过程
列举特征提取、边缘检测算法
- 相机标定介绍
- 理论基础
  - 棋盘格检测
  - 基本符号
  - 图像平面与棋盘格平面之间的单应矩阵
- 计算 A^{-T}A^{-1}
- 计算相机内参矩阵
- 计算相机外参矩阵
  - SVD精调R
  - 优化外参
- 估计镜头的畸变系数
  - 估计畸变的初值
  - 带畸变的代价函数
- 退化配置
- 标定流程总结（包括算法）
1×1 卷积核的作用是什么？
- 附录
  - 矩阵的F范数

齐次坐标、点到直线距离

参考：https://zhuanlan.zhihu.com/p/26307123

直线的一般式：

Ax By C=0 (A、B不同时为0)【适用于所有直线】

两垂直直线的斜率积为-1（如果斜率都存在的话）：

k_1 times k_2 = -1

可以推出点(x, y)到直线的距离公式：

d = frac{|Ax By C|}{sqrt{A^2 B^2}}

给三角形三边求面积

参考：https://baike.baidu.com/item/三角形面积公式/8491990

海伦公式，三角形三边分别为a, b, c，求面积S:

S = sqrt{p(p-a)(p-b)(p-c)} \ p = frac{a b c}{2}

简述SIFT特征点检测、描述和匹配的过程

**SIFT(Scale Invariant Feature Transform)**特征对旋转、尺度缩放、亮度变化等保持不变性，是非常稳定的局部特征。 SIFT的主要思路：a)构造图像的尺度空间表示，b)在尺度空间中搜索图像的极值点，c)由极值点建立特征描述向量，d)用特征描述向量进行相似度匹配。检测：通过图像与DOG算子卷积得到一幅二维图像在不同尺度下的尺度空间表示（即DOG图像）；然后通过每个像素点与其三维领域的临近点进行比较，找出DOG局部极值点作为初步的特征点；然后通过曲线拟合（临近信息插补）得到特征点的精确位置，同时也会舍弃那些不明显关键点和边缘响应；接下来利用特征点邻域像素的梯度分布来确定特征点的方向。每个特征点都包含三个信息(x,y,σ,θ)(x,y,σ,θ)，即位置、尺度和方向。描述：描述子将被用来作为目标匹配的依据，所以应具有较高的独特性以保证匹配率。特征描述大致包含三个步骤：1.旋转主方向，即将坐标轴旋转为关键点方向，以确保旋转不变性。2.然后在特征点对应的高斯图像上统计其16∗1616∗16邻域内的方向梯度，将统计向量作为该点的sift描述子。3.特征向量的常规归一化，进一步去除光照变化等影响。匹配：RANSAC(RANdom SAmple Consensus)方法是当前常用的一种特征点对的匹配算法，在OpenCV中我们可以使用暴力匹配（BruteForceMatcher）和FLANN（Fast Library for Approximate Nearest Neighbors，快速最近邻逼近搜索函数库）实现快速高效匹配。另外由于噪声等其他因素影响，有时次优匹配点可能和最优匹配点非常接近（接近0.80.8），实验证明舍弃这些点效果会更好。参考：Introduction to SIFT (Scale-Invariant Feature Transform)；《数字图像处理MATLAB版（左飞著）》12.1章；《OpenCV3编程入门》第11章

列举特征提取、边缘检测算法

边缘检测：sobel、LoG、Canny

特征提取：harris、SIFT、SURF、ORB

相机标定介绍

张在论文中把相机标定方法分为两类：摄影测量标定法(Photogrammetric calibration) 和 自标定法(Self-calibration)，

摄影测量标定法(Photogrammetric calibration) ：采用已知尺寸的高精度3D标定物进行标定的方法。该类方法标定结果准确，但是需要昂贵的标定物以及复杂的标定设置。属于该类方法的一些工作如下：
- O. Faugeras. Three-Dimensional Computer Vision: a Geometric Viewpoint. MIT Press, 1993.
- R. Y. Tsai. A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf tv cameras and lenses. IEEE Journal of Robotics and Automation, 3(4):323–344, Aug. 1987.
自标定法(Self-calibration) ：不需要任何标定物，只需要在静态场景中移动相机。该类方法标定结果并不总是可靠的，但是标定比较灵活。属于该类方法的一些工作如下：
- S. Bougnoux. From projective to euclidean space under any practical situation, a criticism of self-calibration. In Proceedings of the 6th International Conference on Computer Vision, pages 790–796, Jan. 1998.
- Q.-T. Luong. Matrice Fondamentale et Calibration Visuelle sur l’Environnement-Vers une plus grande autonomie des systemes robotiques. PhD thesis, Universitede Paris-Sud, Centre d’Orsay, Dec. 1992.
- Q.-T. Luong and O. Faugeras. Self-calibration of a moving camera from point correspondences and fundamental matrices. The International Journal of Computer Vision, 22(3):261–289, 1997.
- S. J. Maybank and O. D. Faugeras. A theory of self-calibration of a moving camera. The International Journal of Computer Vision, 8(2):123–152, Aug. 1992.

根据这两种分类原则，张认为自己提出的方法介于二者之间，其采用高精度的2D标定平面（棋盘格），并且只需要移动相机或者移动棋盘格。采用的标定物便宜且易于获取（打印棋盘格在纸张上即可），标定配置灵活（移动棋盘格平面，拍摄各个方向的棋盘格即可），且标定结果鲁棒精度较高（实验表明仿真和真实数据表现都不错）。

其实，除了这两类方法，张还提到了一些其他方法，不太好分类，例如：

vanishing points for orthogonal directions [3, 14], | 这个标定内参效果好像不错，就是vanishing points不好找。
- B. Caprile and V. Torre. Using Vanishing Points for Camera Calibration. The International Journal of Computer Vision, 4(2):127–140, Mar. 1990.
- D. Liebowitz and A. Zisserman. Metric rectiﬁcation for perspective images of planes. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 482–488,Santa Barbara, California, June 1998. IEEE Computer Society.
calibration from pure rotation [11, 21].
- R. Hartley. Self-calibration from multiple views with a rotating camera. In J.-O. Eklundh, editor, Proceedings of the 3rd European Conference on Computer Vision, volume 800-801 of Lecture Notes in Computer Science, pages 471–478, Stockholm, Sweden, May 1994. Springer-Verlag.
- G. Stein. Accurate internal camera calibration using rotation, with analysis of sources of error. In Proc. Fifth International Conference on Computer Vision, pages 230–236, Cambridge, Massachusetts, June 1995.

理论基础

棋盘格检测

在标定之前，需要得到标定板上棋盘格的世界坐标和投影到图像上的像素坐标。因为标定板的尺寸已知，我们将坐标系原点放到标定板的某一点上，所以棋盘格世界坐标是已知的；而棋盘格的像素坐标，是通过棋盘格的角点检测算法检测得到的。

角点检测算法也是一个大领域，且不是张氏标定法的重点，本文不打算展开讨论。

基本符号

相机拍摄已知尺寸的棋盘格，得到的投影点记为 m=[u, ; v]^T ，对应的3D点记为 M=[X, ; Y, ; Z]^T ，采用齐次坐标表示，对应的齐次坐标为 widetilde{m} =[u, ; v, ; 1]^T 和 widetilde{M}=[X, ; Y, ; Z, ; 1]^T ，根据针孔相机模型，我们可以知道他们满足如下关系：

swidetilde{m} = A[R ; t]widetilde{M}

其中，s为标量，表示尺度， [R ; t] 为相机外参，表示从世界坐标系转到相机坐标系下所需要进行的变换，A为相机内参矩阵，有如下形式：

A = begin{bmatrix} alpha & gamma & u_0\ 0 & beta & v_0\ 0 & 0 & 1 end{bmatrix}

其中，alpha, beta 分别是像素坐标系u和v坐标系方向像素焦距， u_0, v_0 是相机主点坐标(即相机光心在图像平面上的投影)，而 gamma 表示像素坐标系u和v之间的倾斜程度(skewness)，在现在的相机中，由于工艺的改善，u和v坐标轴几乎垂直，所以 gamma=0 。

图像平面与棋盘格平面之间的单应矩阵

不失一般性，为了简化推导过程，我们假设标定板在 Z=0 平面上

因此投影关系变为

sbegin{bmatrix} u \ v \ 1 end{bmatrix} = A [r_1 ; r_2 ; r_3 ; t]begin{bmatrix} X \ Y \ Z \ 1 end{bmatrix} = A [r_1 ; r_2 ; t]begin{bmatrix} X \ Y \ 1 end{bmatrix}

为了方便，这里仍然用 M 去表示世界坐标系的点，只不过由于 Z=0 ，这里只关注X和Y，所以 M=[X, ; Y]^T, ; widetilde{M}=[X, ; Y, ; 1]^T ，因此上述关系可以记为

begin{aligned} swidetilde{m} &= Hwidetilde{M} \ H &= A [r_1 ; r_2 ; t] end{aligned} tag{2}

其中，这里的 H 其实就是两个平面之间的单应矩阵。

单应矩阵 H 的求法比较固定，需要4对匹配点即可计算出来，这里就不展开讨论了，假设你已经求出单应矩阵 H ，下面讨论如何根据 H 的结果计算出相机内外参数。(H 的求法，可以看论文的附录A，用的最大似然估计直接进行优化得到，另外4点求解的方法也可以看2015年大疆算法工程师笔试题)

计算 A^{-T}A^{-1}

本节将利用旋转矩阵是个正交矩阵的性质，去计算相机内参。

begin{aligned} r_1^Tr_2 &= 0 \ ||r_1||_2^2 &= ||r_2||_2^2 = 1 end{aligned}

即

begin{aligned} h_1^TA^{-T}A^{-1}h_2 &= 0 \ h_1^TA^{-T}A^{-1}h_1 &= h_2^TA^{-T}A^{-1}h_2 = 1 end{aligned}

将 A^{-T}A^{-1} 记为B：

B = A^{-T}A^{-1} = begin{bmatrix} B_{11} & B_{12} & B_{13} \ B_{12} & B_{22} & B_{23} \ B_{13} & B_{23} & B_{33} end{bmatrix} = =left[begin{array}{ccc} frac{1}{alpha^{2}} & -frac{gamma}{alpha^{2} beta} & frac{v_{0} gamma-u_{0} beta}{alpha^{2} beta} \ -frac{gamma}{alpha^{2} beta} & frac{gamma^{2}}{alpha^{2} beta^{2}} frac{1}{beta^{2}} & -frac{gammaleft(v_{0} gamma-u_{0} betaright)}{alpha^{2} beta^{2}}-frac{v_{0}}{beta^{2}} \ frac{v_{0} gamma-u_{0} beta}{alpha^{2} beta} & -frac{gammaleft(v_{0} gamma-u_{0} betaright)}{alpha^{2} beta^{2}}-frac{v_{0}}{beta^{2}} & frac{left(v_{0} gamma-u_{0} betaright)^{2}}{alpha^{2} beta^{2}} frac{v_{0}^{2}}{beta^{2}} 1 end{array}right]

注意到B是一个对称矩阵，只有6个未知数 ，记为列向量b，即

b = [B_{11}, B_{12}, B_{13}, B_{22}, B_{23}, B_{33}]^T

h_i^TBh_j = v_{ij}^Tb

其中

v_{ij} = [h_{1i}h_{1j}, ; h_{1i}h_{2j} h_{2i}h_{1j}, ; h_{2i}h_{2j}, ; h_{3i}h_{1j} h_{1i}h_{3j}, ; h_{3i}h_{2j} h_{2i}h_{3j}, ; h_{3i}h_{3j}]^T

所以每一张标定板图片带来的两个约束可表示如下：

计算相机内参矩阵

下面我们要通过B求解A了。首先展开B矩阵，有：

B = lambda A^{-T}A^{-1} = lambda left[begin{array}{ccc} frac{1}{alpha^{2}} & -frac{gamma}{alpha^{2} beta} & frac{v_{0} gamma-u_{0} beta}{alpha^{2} beta} \ -frac{gamma}{alpha^{2} beta} & frac{gamma^{2}}{alpha^{2} beta^{2}} frac{1}{beta^{2}} & -frac{gammaleft(v_{0} gamma-u_{0} betaright)}{alpha^{2} beta^{2}}-frac{v_{0}}{beta^{2}} \ frac{v_{0} gamma-u_{0} beta}{alpha^{2} beta} & -frac{gammaleft(v_{0} gamma-u_{0} betaright)}{alpha^{2} beta^{2}}-frac{v_{0}}{beta^{2}} & frac{left(v_{0} gamma-u_{0} betaright)^{2}}{alpha^{2} beta^{2}} frac{v_{0}^{2}}{beta^{2}} 1 end{array}right]

这里有6个已知数 [B_{11}, B_{12}, B_{13}, B_{22}, B_{23}, B_{33}] ，也有6个未知数 [lambda, alpha, beta, gamma, u_0, v_0] ，可以按下列顺序依次求解出所有未知数：

begin{aligned} v_0 &= (B_{12}B_{13}-B_{11}B_{23})/(B_{11}B_{22}-B_{12}^2) \ lambda &= B_{33}-[B_{13}^2 v_0(B_{12}B_{13}-B_{11}B_{23})]/B_{11} \ alpha &= sqrt{lambda / B_{11}} \ beta &= sqrt{lambda B_{11} / (B_{11}B_{22}-B_{12}^2)} \ gamma &= -B_{12} alpha^2 beta / lambda \ u_0 &= gamma v_0 / beta - B_{13} alpha^2 / lambda end{aligned}

上面公式是论文里给出来的，我只验证了前两个是对的，后面可能有笔误，建议自己动手推一下。

至此，内参矩阵 A 的所有未知参数已计算出来。下面利用 A 计算相机外参。

计算相机外参矩阵

当相机内参矩阵 A 已知时，由前面公式 (2) 可以求解出相机外参：

r_1 = lambda A^{-1}h_1 \ r_2 = lambda A^{-1}h_2 \ r_3 = r_1 times r_2 \ t = lambda A^{-1}h_3

这里的 lambda = 1/||A^{-1}h_1|| = 1/||A^{-1}h_2|| 。r_1, r_2, r_3 可以恢复出旋转矩阵 R ，但由于数据带有噪声，恢复的 R = [r_1, r_2, r_3] 通常不满足旋转矩阵的性质(即R为正交矩阵，且行列式为1)。所以在恢复出R后，还需要通过 SVD 分解的方式在 R 附近，去找到一个满足旋转矩阵性质的 R 。

SVD精调R

通过SVD精调R的方法很常用，你会在许多算法的推导里看到，建议弄懂。

为了不混淆，这里将上面求解出来的，不一定满足旋转矩阵性质的解记为 Q ，而将我们期待得到的，符合性质的旋转矩阵记为 R 。那么寻找这么一个矩阵 R 的问题是一个最优化问题，代价函数为：

underset{R}{min} ||R-Q||_F^2, quad subject ; to ; R^TR=I. tag{15}

因为

begin{aligned} ||R-Q||_F^2 &= trace((R-Q)^T(R-Q)) \ &= trace(R^TR Q^TQ - 2R^TQ) \ &= 3 trace(Q^TQ) - 2trace(R^TQ) end{aligned}

注意：注意：这里虽然看上去和ICP那里很像，但是推导依据不一样，ICP是根据正定矩阵迹的性质推的，而这里Q不是正定矩阵。

至此，我们通过SVD找到了更好的旋转矩阵R(符合旋转矩阵的性质)。

优化外参

前面已经计算出了相机内外参，但是一方面由于数据存在噪声，另一方面在计算R时通过SVD找了个最近的解，所以得到的相机内外参并不一定能使得相机匹配点的投影误差最小，可以通过最优化的方法进一步调整内外参。

假设拍摄了n张标定板的图片，每张图片上棋盘格的角点个数有m个，则代价函数如下：

sum_{i=1}^{n}sum_{j=1}^{m} ||m_{ij} - hat{m}(A, R_i, t_i, M_j)||^2

其中，hat{m}(A, R_i, t_i, M_j) 是 M_j 通过公式 (2) 投影到图像上的投影点。

张的方法中，将R用旋转向量表示，来进行优化。我个人觉得旋转向量优化时还是会出问题，因为存在周期，360°等于0°。用李代数可能更好。

具体优化过程略，可以参考常见的最优化方法，如最速下降法、高斯-牛顿法、LM算法。

估计镜头的畸变系数

实际情况中相机总是要配备镜头，而镜头都是带有畸变的，前面的讨论都没有涉及到镜头畸变，所以算出来的内外参，哪怕优化后也是有误差的。

估计畸变系数的方法，论文里没有明确说，个人理解大概能分成三种：

反复迭代的方法。先不考虑畸变，计算出内外参，再通过内外参求解畸变系数，再考虑畸变系数，计算出内外参，再通过内外参求解畸变系数，再考虑畸变系数…一直反复迭代到收敛。张的论文中提到这种方法在实验过程中收敛较慢。
最优化方法。在代价函数中添加畸变系数项，畸变系数初始值设置为0，把前面算出来的内外参作为初始值丢进去，一起优化。
结合上述两种方法。在代价函数中添加畸变系数项，先不考虑畸变，计算出内外参，再通过内外参求解畸变系数，此时的畸变系数和内外参作为优化问题的初始值，优化。该方法与第二种方法的区别就是，畸变系数的初始值不是0，而是先大概计算了一下初始值。

这里直接讨论第三种方法。

估计畸变的初值

首先估计畸变的初值。论文中说镜头畸变包括径向畸变和轴向畸变以及一些其他畸变，但起主要作用的，主要是径向畸变的前两项 k_1, k_2 ，考虑更多的畸变参数并不一定能得到更好结果，反而可能因为数值稳定性问题使得去畸变效果变差。

论文仅仅考虑径向畸变的前两项，畸变公式为：

hat{x} = x x[k_1(x^2 y^2) k_2(x^2 y^2)^2] \ hat{y} = y y[k_1(x^2 y^2) k_2(x^2 y^2)^2]

这里， k_1, k_2 为径向畸变的前两项系数， (x, y) 为无畸变的相机成像平面坐标(即归一化平面坐标，不是像素坐标)，(hat{x}, hat{y}) 是镜头畸变后的相机成像平面坐标。

根据归一化平面坐标与像素坐标的变换关系：

hat{u} = u_0 alphahat{x} gammahat{y} \ hat{v} = v_0 betahat{y}

畸变公式在像素坐标系下表示为：

hat{u} = u (u-u_0)[k_1(x^2 y^2) k_2(x^2 y^2)^2] \ hat{v} = v (v-v_0)[k_1(x^2 y^2) k_2(x^2 y^2)^2]

为了计算未知数 k_1, k_2 ，将畸变公式改写成 Ax=b 的形式：

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} hat{u}-u \ hat{v}-v end{bmatrix}

其中， (u, v) 是将标定板坐标通过前面估计的内外参投影得到， (hat{u},hat{v}) 是对标定板图像执行棋盘格检测得到的像素坐标。

对于每一个像素点，都能得到上述的两个约束方程，假设我们拍摄了n张标定板的图像，每张图像上有m个点，那么我们可以得到2nm个约束方程，我们将系数矩阵记为D，有

Dk=d

其中，D为 2nm times 2 大小的系数矩阵， k = [k_1, ; k_2]^T ， d 为 2nm times 1 大小的列向量。

求解该方程组的方法有很多，具体可以看我写的另一篇文章矩阵分解与线性方程组。张的论文里是通过正规方程的方式求解的：

k = (D^TD)^{-1}D^Td

求解该方程组后，我们就可以得到畸变系数 k_1, k_2 ，将这作为初值，与前面计算得到的内外参一起进行最优化即可。

带畸变的代价函数

同前面的优化部分一样，这里只给出代价函数，具体优化过程略，可以参考常见的最优化方法，如最速下降法、高斯-牛顿法、LM算法。

假设拍摄了n张标定板的图片，每张图片上棋盘格的角点个数有m个，则代价函数如下：

sum_{i=1}^{n}sum_{j=1}^{m} ||m_{ij} - hat{m}(A, k_1, k_2, R_i, t_i, M_j)||^2

其中，hat{m}(A, k_1, k_2, R_i, t_i, M_j) 是 M_j 通过公式 (2) 投影到图像上的投影点。

张的方法中，将R用旋转向量表示，来进行优化。我个人觉得旋转向量优化时还是会出问题，因为存在周期，360°等于0°。用李代数可能更好。

其中对于畸变系数 k_1, k_2 的初值，可以给0，也可以给前一小节中估计出来的值。当然，给的初值越好，优化收敛速度越快，且越不容易收敛到错误的极小值中去。

退化配置

数学推导请看原论文，以后有时间再补吧。

一句话概括：若有标定板之间平行，则只有一个标定板提供了约束，因为其他标定板的单应矩阵，都能由这个标定板单应矩阵的列向量线性组合得到。所以若想得到好的标定结果，则标定板之间不应该平行。

个人看法：虽然说多余的平行标定板没有用，但这是对于无畸变图像来说的，对于有畸变图像来说，在标定相机内外参初值时用处可能较小，但在估计畸变时，是有用的。另外，越靠近图像边缘，畸变越明显，如果因为某些原因仅能拍摄少数标定板图像，则建议拍摄不平行的，位于画面四周的标定板（画面正中心拍摄一两张应该就可以了）。

标定流程总结（包括算法）

标定板采用的是棋盘格，标定板棋盘格尺寸已知。标定流程总结如下：

拍摄多张不同方向的标定板图像，并进行棋盘格角点检测
对于每张图像，根据检测到的角点信息和棋盘格尺寸信息，都可以计算出一个对应的单应矩阵 H
利用单应矩阵 H 计算 B 矩阵 B = A^{-T}A^{-1}
求解出 B 矩阵后，进一步求解得到内参 A
求解出 A 矩阵后，进一步求解得到外参 R 和 t
根据 A, R, t 估计出相机畸变系数(初值)
构建带畸变系数的最优化问题，同时优化 A, R, t 直到收敛

1×1 卷积核的作用是什么？

修正线性激活（ReLU）
实现跨通道的交互和信息的整合
进⾏卷积核通道数的降维和升维

实现多个特征映射（feature map）的线性组合，实现通道个数的变换
对特征图像进⾏⽐例缩放
减少参数量

附录

矩阵的F范数

和向量的2范数类似，矩阵也有2范数，只不过矩阵中叫做 m_2 范数，或者叫F范数。对于 m times n 大小的矩阵 A ，其F范数定义如下：

||A||_F = (sum_{i=1}^{m} sum_{j=1}^{n} |a_{ij}|^2)^{1/2} = (trace(A^HA))^{1/2}

其中 A^H 表示是 A 的 共轭转置 矩阵，若 A 是实矩阵，则 A^H = A^T 。

三维视觉面试，复习笔记：

https://github.com/Liber-coder/CV_Notes

图像处理100问：

https://github.com/gzr2017/ImageProcessing100Wen
Zhang, Zhengyou. “A flexible new technique for camera calibration.” IEEE Transactions on pattern analysis and machine intelligence 22.11 (2000): 1330-1334.
张正友标定算法原理详解 | 不错的资料，借鉴了部分

camera self

0 人点赞