【草履虫都能学会】02 计算机视觉数学基础

前言

  01篇介绍了相关的物理基础和几何变换基础，本篇接着介绍计算机视觉需要掌握的数学基础。

一、线性时不变系统

摄影测量学涉及一系列图像处理技术，例如平滑、去噪等。把这些相关图像处理技术纳入到“线性时不变系统”这个统一的数学体系下。

什么是线性时不变系统？ 顾名思义，线性+时间不变。线性系统指的是满足线性叠加原理的输入输出系统，也就是满足乘法交换律和结合律。下式就是满足线性系统的一个例子：
$$L(a{x}_1+b{x}_2)=aL(x_1)+bL(x_2)$$而时不变系统，指的是输入和输出之间的作用不受时间起始点的影响（随着时间不会发生变化，例如太阳无论如何都是东升西落）。同时满足两者的称为线性时不变系统。

一个线性时不变系统一定能够用卷积表达，反之也成立，即能用卷积表达的系统也一定是线性时不变系统。这样卷积就有了线性时不变系统的所有特性，例如结合律：
$$h\ast (g\ast f)=(h\ast g)\ast f$$

二、相关应用

介绍了卷积的相关实际应用，帮助读者更好的理解卷积操作以及线性时不变系统在图像处理中的应用。此处不需要特别熟练掌握，但是要知道讲的是什么事，更多是先培养相关的基础思维。

1. 边缘提取

边缘提取事实上就是将一个对边敏感的模板作用于原始图像上，例如Canny算子模版，Sobel算子模版等。把模版在图像上滑动并计算，从而进行图像处理。算子自身的特性决定了计算的方式，同样决定了对于边缘的敏感性。
例如上图为Sobel算子，设计了3x3的卷积核，通过与图像进行卷积计算，能够得到图像的水平梯度和垂直梯度。梯度实际就是原始图像的极值表达，对于边缘部分，灰度会有急速的变化，Sobel算子能够检测出这种变化，从而检测出边缘。梯度计算的表达式如下：$$J(\mathbf{x})=\nabla I(\mathbf{x})=\left(\frac{\partial I}{\partial x},\frac{\partial I}{\partial y}\right)$$式中的$I(x)$可以理解为待处理的原始图像，$J(x)$代表计算得到的梯度图。不过注意在离散图像中，一般用差分来代替微分。另外，通常需要在计算梯度之前先进行高斯滤波去噪，用一个高斯卷积核$G_{\sigma }$作用于图像$I(x)$，如下：
$$J_{\sigma }(\mathbf{x})=\nabla (G_{\sigma }\ast I(\mathbf{x}))=(\nabla G_{\sigma })(\mathbf{x})\ast I(\mathbf{x})$$

有时候，会进一步计算梯度图的极值点，也就是计算梯度的梯度，从而得到更特殊的边缘。这种求取梯度的梯度操作即为LapLacian，其计算公式如下：
$$\nabla J_{\sigma }(\mathbf{x})=\nabla ((\nabla G_{\sigma })(\mathbf{x})\ast I(\mathbf{x}))=({\nabla }^2{G}_{\sigma })(\mathbf{x})\ast I(\mathbf{x})$$式子中的卷积核${\nabla }^2{G}_{\sigma }$又被称为LapLacian of Gaussian(LoG)，其实就是对高斯卷积核做了二阶梯度，从而结合高斯平滑和边缘检测。公式表达如下：
$${\nabla }^2{G}_{\sigma }(\mathbf{x})=\left(\frac{{\partial }^2{G}_{\sigma }}{\partial x^2},\frac{{\partial }^2{G}_{\sigma }}{\partial y^2}\right)$$LoG算子的应用十分广泛，它是最低阶的各向同性算子（各向同性即旋转不变性，即可以在存在旋转的情况下完成边缘检测）。

涉及到的公式不用记忆，只是为了帮助读者更好的理解，请不要有数学恐惧。

2. 图像金字塔

什么是图像金字塔？可以理解为一幅图像的不同尺度的组合，尺度越大的图像在金字塔的越底层。例如下图中，最底层的图像尺度最大，相应的图像细节越多，图像分辨率越高，越清晰。随着尺度的减小，图像也就越来越模糊，最顶层的图像最模糊，细节丢失最多。

如何获取图像金字塔？利用高斯卷积核。高斯滤波可以对图像进行模糊平滑，每进行一次平滑就会让图像变模糊一个层次。例如先对最底层的原始图像进行高斯滤波，然后再每隔2个像素进行采样，用采样得到的像素生成高一层的金字塔，以此类推。这样得到的图像每隔一层长和宽会减少一半，从而得到图像金字塔。

读者可能会疑惑，为什么不直接采样，还要高斯滤波？因为直接采样会存在锯齿效应，因此先进行高斯平滑，让图像变模糊一点，这样采样的时候就不会有明显的锯齿效应，当然图像随层数增加也会变得越来越模糊。

建立图像金字塔的作用。金字塔在图像匹配中有重要作用，借助高斯金字塔，可以实现由粗到精的
搜索策略，加快图像匹配速度（系列后续会详细介绍）。能在最佳分辨率上寻找和定位感兴趣的物体或特征。能够满足超大图像的实时浏览，例如读取遥感图像需要较大的内存容量，而利用高斯金字塔，就可以实现大尺度图像的压缩，节省内存空间。

3. 图像匹配

寻找两张图像中的相同物体或相同特征。在遥感中通常是进行两种相同地物的匹配。在计算机视觉三维重建中，图像匹配可以理解为两张图像中同名特征点的相互匹配。具体会在后续文章中介绍，此处先略过。

三、最小二乘平差

最小二乘平差可真是重量级内容，用到的地方超级多，务必仔细阅读学习，做到熟练掌握。公式也必须了解熟悉，推导过程也得知道一二，最好在学习之后能够自己推导出来

为什么引入最小二乘平差？因为有个东西叫做“最优化问题”，打过数学建模的同学可能比较了解。最优化问题几乎出现在所有的科学和工程领域，是一种将观测值误差和模型误差进行最佳分配的策略。伟大数学家和大地测量学家高斯在19岁就提出了最小二乘平差，在测绘领域中沿用至今（我只能说高斯YYDS）。最小二乘平差就是为了解决最优化问题提出的一种方法。

1. 前提条件

使用最小二乘平差有两个前提条件：

1. 观测值多于未知数。 也就是需要进行多余观测，完成的是超定方程组的求解。
2. 观测值误差是偶然误差（随机误差），且误差服从高斯分布 。注意区分偶然误差和系统误差，偶然顾名思义就是偶然出现的随机出现的，不可控的误差。在摄影测量中，一般将偶然误差视作服从高斯分布，也叫做随机误差。

2. 什么是二乘？

这里要知道范式的定义：$\Vert x{\Vert }^p$就是p范式，展开表示为：
$$\begin{array}{rl}\Vert x{\Vert }^p& =\left(\sum _{i=1}^n|x_i{|}^p\right)=|x_1{|}^p+|x_2{|}^p+\cdots +|x_n{|}^p\end{array}$$那么二次范式$L_2$（p=2）就是二乘的意思，最小二乘就是平方和最小的意思。 $min\Vert x{\Vert }^2$

至于为什么选择二乘，不选一乘，三乘？因为$L_2$下的最优化问题可以转化为为解线性方程组，解空间是凸的（vertex），并有一个全局最优解。但其他范式的解空间是非凸的（non-vertex），不但需要非线性迭代，同时可能陷入局部最优解。

3. 最小二乘平差定义

设x为n维未知变量，$F(x)=[f_1(x),...,f_m(x){]}^T$是关于x的函数，我们需要找到一组最佳的x，使得下式成立$$\underset{x}{\min }\frac{1}{2}\Vert F(x){\Vert }^2$$如果$F(x)$是线性的，则可以写为$F(x)=Ax+b$，则$\Vert F(x){\Vert }^2$可以看做复合二次函数，要求其最小，需要对$\Vert F(x){\Vert }^2$求导找极值。求导得导数为0位置$$A^T(Ax+b)=0$$也就是$$A^{T}Ax=-A^{T}b$$求解上式线性方程组即可解值。

4. 不等权的最小二乘平差

之前假设所有观测值x都是等精度的，其函数$f_i(x)$也是等精度的。但是某些情况下最优化问题可能有精度不同的观测值，则需要更普遍的目标函数，进行加权可得$$\underset{x}{\min }\frac{1}{2}P\Vert F(x){\Vert }^2$$P是权矩阵，对角线矩阵（通常假定观测值之间相互独立），求导得$$A^{T}PAx=-A^{T}Pb$$

5. 非线性最小二乘平差（高斯-牛顿法）

前两种都假设观测值函数是线性的，但更多情况下其实是非线性的，例如透视变换。此时需要将非线性的$F(x)$函数转换为线性，然后再进行平差求解。

对$F(x)$泰勒级数展开至一次项，即$F(x+\Delta x)\approx F(x)+J(x)\Delta x$，$J(x)$是$F(x)$的一阶导数，目标函数为$$\underset{x}{\min }\frac{1}{2}\Vert F(x)+J(x)\Delta x{\Vert }^2$$解为$$J(x{)}^{T}J(x)\Delta x=-J(x{)}^{T}F(x)$$但是由于泰勒展开截断了二次以上的项，所以一次近似解算不够，一般需要迭代：用解得的$\Delta x$不断修正$x$，即用$x^{k+1}=x^k+\Delta x^k$计算新的$J(x)$和$F(x)$，然后再次带入上式求解，直到满足限定的条件（迭代结束条件），例如满足迭代前后解相差不大$$\frac{\Delta x_{k+1}^i}{\Delta x_k^i}\le T$$其中上标$i$代表$\Delta x$的任意一个分量，一般$T$取$10^{-5}$。这种非线性的最小二乘迭代解法，在数学上叫做高斯-牛顿法。

6. Hessian矩阵

Hessian矩阵用$H$表示，定义为：$$H(x)=J(x{)}^{T}J(x)$$该矩阵$H(x)$称为信息矩阵，Hessian矩阵，设计矩阵，法方程系数矩阵等。而$g=J(x{)}^{T}F(x)$被称为法方程常数项，是目标函数$\underset{x}{\min }\frac{1}{2}\Vert F(x){\Vert }^2$的一阶导数（梯度）。

7. LM算法

正常情况下，可以对高斯-牛顿法中$J(x{)}^{T}J(x)\Delta x=-J(x{)}^{T}F(x)$求逆解$\Delta x$，但是如果$H(x)$的对角线元素为0或接近0，是奇异矩阵或接近奇异矩阵（$|H(x)|$=0，此时H不可逆），那么解方程就不能求逆或求逆变成了病态问题。

于是需要引入正则化项，对H矩阵的对角线作进一步约束。此时目标函数变为$$\underset{x}{\min }\frac{1}{2}\Vert F(x)+J(x)\Delta x{\Vert }^2+\mu \Vert D(x)\Delta x{\Vert }^2$$对应的解为$$[J(x{)}^{T}J(x)+\mu D(x{)}^{T}D(x)]\Delta x=-J(x{)}^{T}F(x)$$其中$D(x)$是$H$矩阵对角线元素的根方差，$\mu$控制$\Delta x$的步长，$\mu$越小，则步长越大，$\mu =0$即传统高斯最小二乘。这种方法被称作 Levenberg-Marquardt (LM)算法。

至此，最小二乘平差的体系介绍完毕。这是针对$L_2$范式下的最优化所发展起来的一套完善的数学方法，广泛应用于摄影测量学。请读者务必仔细揣摩，熟练掌握。

四、粗差和系统误差的处理

1. 粗差处理——RANSAC

最小二乘并不是抗差算法。抗差指抵抗粗差的能力，而粗差是明显偏离数学模型的观测。最小二乘假定观测值符合高斯分布；在现实条件下，该条件往往无法满足。以最直观的直线拟合为例，若观测值中只混进一个粗差，采用最小二乘平差将得到明显错误的结论。对粗差的恰当处理，是保证最终平差结果正确性的先决条件。

这里介绍一种抗差算法：RANSAC，又称为“随机采样一致性检测方法”。

假设拟合一条直线，可以由两个点确定，模型参数$m=2$。随机抽取两个点（最小样本点），若这两个点都是内点（符合模型的观测值），则这两个点可以确定正确的直线$l$，把其余多余观测点带入该直线，如果是内点，则符合$x\cdot l=0$（点在线上），而粗差则不满足该方程，可以被探测出来。这种探测粗差并同时算出正确模型参数的方法称为RANSAC。其核心在于确定内点。

令观测值个数为$n$，必要观测数（确定模型需要的最小观测数）为$m$，其中$n\gg m$。设内点的比例为$w$，每次随机抽取$m$个点，抽取$k$次的情况下，得到$m$个点恰好都是内点的概率$p$ $$p=1-(1-w^m{)}^k$$有理论研究在粗差高达75%的情况下，迭代一定次数后依然能保证99%的正确率，所以RANSAC是一个良好的抗差算法。

算法步骤：

1. 以直线拟合为例，抽样$k$次（给$k$赋值开始循环）。
2. 每次随机抽$m$个点，计算一次直线参数得到直线方程.
3. 把剩余的点代入直线方程，如果值为0，则判断为内点，并计算内点数量。不为0则为粗差。
4. 循环结束后，内点数目最多的那个直线方程，就是最佳的方程，也就是最终的模型。

缺点： 需要采样，导致只能处理低维参数空间中的粗差，也不能处理时间序列等动态模型中的粗差。

优点：优点当然显而易见。 原理简单、易嵌于入，在摄影测量和计算机视觉中应用广泛。

2. 系统误差

系统误差，也叫模型误差。顾名思义，模型误差就是模型自身不能代表物理真实而呈现的系统性偏差。例如，在大质量恒星附近，牛顿光学（认为光沿直线传播）将产生可观的偏差，我们就需要相对论这个更加严格的模型。

第二个例子是相机检校。在系列后续会细讲，简单说就是由于光学镜头的系统误差导致像片发生畸变，于是需要利用多项式模型进行畸变的修正。

第三个例子是卫星摄影测量的成像模型。卫星摄影测量利用的是线阵CCD推扫成像，不同与普通相机的面阵CCD一次性成像。由于卫星成像的特殊性，我们使用RFM（有理多项式模型）代替卫星成像模型。但这种几何上的强行拟合势必引入一些模型误差；此时，通常在像方坐标中加入一个仿射变换，以修正 RFM 模型，得到更好的模型定位精度。

由上述三个例子可以看出，模型误差处理的两个基本原则为：

1. 建立更严密的物理模型
2. 物理模型较复杂时，采用几何拟合模型，并可以再引入其他几何模型（几何模型之间一般线性叠加）

小结

关于计算机视觉的高等数学基础，需要有一定的微积分基础（知道什么是微分，以及能够求解偏导数，求解全微分）和一定的线性代数基础（知道什么是矩阵，了解矩阵的相关运算）。

另外会涉及少许图像处理的内容，例如卷积操作，滤波和去噪等，读者如果没有基础的话可以先去补补。