图像算子概述

用机器视觉系统分析未知场景时，计算机并不预先知道图像中物体的尺度。我们需要同时考虑图像在多尺度下的描述，获知感兴趣物体的最佳尺度。如果不同的尺度下都有同样的关键点，那么在不同的尺度的输入图像下就都可以检测出来关键点匹配，也就是尺度不变性。

尺度空间中各尺度图像的模糊程度逐渐变大，能够模拟人在距离目标由近到远时目标在视网膜上的形成过程。
尺度越大图像越模糊。

尺度空间表达与金字塔多分辨率表达

高斯核是唯一可以产生多尺度空间的核（Scale-space theory: A basic tool for analysing structures at different scales）。一个图像的尺度空间$L(x,y,σ)$ ,定义为原始图像$I(x,y)$与一个可变尺度的2维高斯函数$G(x,y,σ)$卷积运算。

二维空间高斯函数：

$$G(x_{i},y_{i},\sigma)=\frac{1}{2\pi\sigma}\exp\left( -\frac{(x-x_{i})^2+(y-y_{i})^2}{2\sigma^2} \right)$$

尺度空间：

$$L(x,y,\sigma)=G(x,y,\sigma)\ast I(x,y,\sigma)$$

尺度是自然客观存在的，不是主观创造的。高斯卷积只是表现尺度空间的一种形式。
二维空间高斯函数是等高线从中心成正太分布的同心圆：

S分布不为零的点组成卷积阵与原始图像做变换，即每个像素值是周围相邻像素值的高斯平均。一个5*5的高斯模版如下所示：

高斯模版是圆对称的，且卷积的结果使原始像素值有最大的权重，距离中心越远的相邻像素值权重也越小。
在实际应用中，在计算高斯函数的离散近似时，在大概3σ距离之外的像素都可以看作不起作用，这些像素的计算也就可以忽略。所以，通常程序只计算(6σ+1)*(6σ+1)就可以保证相关像素影响。

高斯模糊另一个很厉害的性质就是线性可分：使用二维矩阵变换的高斯模糊可以通过在水平和竖直方向各进行一维高斯矩阵变换相加得到。

$O(N^2\ast m\ast n)$次乘法就缩减成了$O(N\ast m\ast n)+O(N\ast m\ast n)$次乘法。（N为高斯核大小，m,n为二维图像高和宽）。

金字塔多分辨率

金字塔是早期图像多尺度的表示形式。图像金字塔化一般包括两个步骤：使用低通滤波器平滑图像；对平滑图像进行降采样（通常是水平，竖直方向1/2），从而得到一系列尺寸缩小的图像。

上图中（a）是对原始信号进行低通滤波，（b）是降采样得到的信号。
而对于二维图像，一个传统的金字塔中，每一层图像由上一层分辨率的长、宽各一半，也就是四分之一的像素组成：

多尺度和多分辨率

尺度空间表达和金字塔多分辨率表达之间最大的不同是：

• 尺度空间表达是由不同高斯核平滑卷积得到，在所有尺度上有相同的分辨率；
• 而金字塔多分辨率表达每层分辨率减少固定比率。
  所以，金字塔多分辨率生成较快，且占用存储空间少；而多尺度表达随着尺度参数的增加冗余信息也变多。
  多尺度表达的优点在于图像的局部特征可以用简单的形式在不同尺度上描述；而金字塔表达没有理论基础，难以分析图像局部特征。

---

Moravec算子

1977年，Moravec提出了兴趣点（Points of Interests)的概念，并应用于解决Stanford Cart的导航问题。1981年， Moravec在International Joint Conference on Artificial Intelligence发表了篇题为：Obstacle Avoidance and Navigation in the Real World by a Seeing Robot Rover的文章，并将其应用与立体匹配。

Moravec的原理如果有一句话来说就是：通过滑动二值矩形窗口寻找灰度变化的局部最大值。具体来说主要包括四个过程：

滑动窗口计算灰度变化（Calculate intensity variantion from shifting windows）

滑动窗口在现有的技术中已经有了很多应用，如模板匹配、目标检测（hog特征的行人检测）等。在Moravec算子中，一般窗口的大小取3×3、5×5、7×7等等，但是随着窗口的增大，计算量也就越大。Moravec算子通过对窗口的水平、垂直和对角八个方向进行移动（Horizontally、Vertically and four diagonals），计算原窗口与滑动窗口差的平方和来得到灰度的变化。我们进一步通过下图一个3×3的滑窗来进行说明：

图中，红色框表示的是原始框，而蓝色框表示向右上的滑动框，白色框表示前景255，黑色框表示背景0。那么原始框和滑动框的灰度变化通过对应位置差的平方和来表示，也即通过下式来计算：

$$\sum_{i=1}^{9}{(A_{i}-B_{2})^2}$$

同样，根据上式计算另外七个方向滑动框的灰度变化（水平向左、水平向右、垂直向上、垂直向下以及四个对角）。至此，我们就计算完成了8个方向的灰度变化，我们称此操作位Moravec operator（Moravec算子）。

构造角点性映射图（Construct the cornerness map）

在构造角点映射图之前，我们先来分析下，通过上式的我们可以得到角点吗？或者凭什么通过计算两个框对应位置的差的平方和就可以检测到角点？问题问得好，我们来看下面的图：

上面四张图上的四个红色的框表示我们正在处理的窗，第一幅图中的窗在表示在目标内部或者是背景上，该区域灰度分布均与，通过对其在8个方向上灰度，灰度变化很小；第二幅图中的窗跨在图像的边缘处，当垂直于边缘方向滑动窗口时将会导致灰度的很大变化，而沿着边缘滑动窗时，灰度变化较小；第三幅图中的窗在角点处，不管往哪个方向滑动窗口，都会导致灰度的很大变化；而第四幅图中的框内是一个离散点，滑动窗向任意方向滑动也会导致灰度的很大变化。

因此，通过上面的描述和分析，我们可以知道，Moravec算子可以作为一种角点性的度量，这种度量是通过求8个方向的滑窗来的最小值来表示。用公式表示如下：

$$C(x,y)=min(V_{u,v}(x,y))$$

我们通过下图来描述角点映射图的构造：

上图中的是通过Moravec算子计算得到的，其中1表示$1\times 255^2$，2表示$2\times 255^2$。通过上图可以知道：

1. 角点位于局部最大值处，我们可以应用非极大值抑制找到局部最大值（non-maximal suppression）。
2. 离散点（噪声点）与角点有相同的角点性（cornerness），因此Moravec算子对噪声敏感，但是通过增大滑窗的大小可以对噪声起到一定的抑制作用，可同时增加了计算量。另一方面，可以通过设定一个阈值T来对cornerness map进行二值化，小于阈值T的cornerness map设置为0，从而对离散点的局部最大值进行抑制。
   阈值的选择引用原文的话：

   Choosing this threshold is difficult as it must be set high enough to avoid these false corners(isolated pixel), but low enough to retain as many true corners as possible.

3. Moravec算子不能应用与图像边界的一定区域（标记为X的区域），对于这部分区域，一般直接忽略，在cornerness map中这些区域对应的值置0。

Moravec算子对角点的检测效果还不错，但是对于对角线上的角点容易出现误检。

总结

Moravec算子作为第一个广泛应用的角点检测算法，开创了角点检测的新纪念，后续的很多角点检测算子都是在其基础上通过扩展得到的。但是Moravec算子也存在诸如方向各异性、噪声敏感、对旋转不具备不变形（角点不具备repeatability）、滑动窗内的各个像素权重同质性（中心像素权重大，离中心越远，权重越小）的缺点，有待改进。
三篇参考：
https://blog.csdn.net/kezunhai/article/details/11176065
https://code.google.com/archive/p/jfeaturelib/
https://www.ri.cmu.edu/pub_files/pub4/moravec_hans_1980_1/moravec_hans_1980_1.pdf

---

Harris算子

在前面的博文中，介绍了Moravec算子，并对moravec算子的不足也进行了简单的描述，Harris算子针对Moravec算子的不足，提出了以下的改进……

Moravec算子各向异性响应（Anisotropic Response of Operator）

Moravec算子仅仅在8个方向（水平、垂直和四个对角方向）计算灰度变化，为了对其扩展，有必要设计一个可以在任何方向对灰度变化进行测度的函数。1988年，Harris和Stephen通过对Moravec算子进行展开，推导得到了Plessey算子，也即Harris算子。

我们先来看看与Harris相关的背景知识。通常，Prewitt算子被用来对图像的梯度进行近似。然而，在实际应用中，一阶梯度通过下图中的公式来进行近似：

对Morevec算子进行分析可以得到：Two Morevec windows中对应像素差的和可以作为图像梯度的合理近似。我们再来看下图：

通过对上图的分析，我们有可以进一步得到：morevec算子中的灰度变化可以采用图像梯度进行近似。灰度的变化可以表示为图像梯度的函数，公式表示如下：

$$V_{u,v}(x,y)=\sum_{\forall i \in mask(x,y)}^{}{\left( u\frac{\partial I_{i}}{\partial x} +v\frac{\partial I_{i}}{\partial y}\right)^2}$$

其中，（u,v)表示滑动，x方向为（1,0），y方向为（0,1），微分的计算如上图所示。

到这里，大家非常明了：上式可以对moravec算子中的灰度变化计算进行精确的逼近。但是又与Moravec算子中灰度变化不同的是通过合理的选择（u，v）可以对任何方向的灰度变化进行测度。

噪声响应（Noise Response）

在Moravec算子中，滑动窗采用的是方形的（square window），方形窗使得不同方向上的中心像素与边界像素的欧式距离是变化的。为了克服这个问题，Harris&Stephen提出只需将方向窗改成圆窗（circle window）。同时，窗中的每个像素是同等地位的，理论上应该是离中心越近的权重越大，而离中心越远，权重越小，因此我们加入高斯权重。因此，灰度变化的新测度方式可以通过下图来表示：

通过公式表示如下：

$$V_{u,v}(x,y)=\sum_{\forall i \in mask(x,y)}^{}{W_{i}\left( u\frac{\partial I_{i}}{\partial x} +v\frac{\partial I_{i}}{\partial y}\right)^2}$$

其中，wi表示位置i处的高斯权重。

边缘的强响应（Large Response of Edge）

因为Moravec算子在边缘处很容易出现误检，Harris&Stephen通过考虑不同方向的灰度度量形成新的角度性测度（cornerness measure）。接着，我们对上面的式子进行变换，如下式：

$$\begin{equation} \begin{split} V_{u,v}(x,y)&=\sum_{\forall i\in mask(x,y)}^{}{\omega_{i}\left( u\frac{\partial I}{\partial x} +v\frac{\partial I}{\partial y}\right)^2}\\ &=\sum_{\forall i\in mask(x,y)}^{}{\omega_{i}\left( u^2\frac{\partial I^2}{\partial x} +2uv\frac{\partial I}{\partial x}\frac{\partial I}{\partial y}+v^2\frac{\partial I^2}{\partial y}\right)}\\&=Au^2+2Cuv+Bv^2 \end{split} \end{equation}\\ where A=\left( \frac{\partial I}{\partial x} \right)^2\otimes\omega_{i}，B=\left( \frac{\partial I}{\partial x} \right)^2\otimes\omega_{i}，C=\left(\frac{\partial I}{\partial x}\frac{\partial I}{\partial y}\right)\otimes\omega_{i},$$

Harris&Stephen同时也注意到，上式可以写成：

$$\begin{equation} \begin{split} V_{u,v}(x,y)&=Au^2+2Cuv+Bv^2 \\&=\left[\begin{matrix}u&v\end{matrix}\right]M\left[\begin{matrix}u\\v\end{matrix}\right]\end{split} \end{equation}\\ where：M=\left[\begin{matrix}A&C\\C&B\end{matrix}\right]$$

对上面的矩阵M，其特征值与图像表面的主曲率是成正比的，并且形成了对M的旋转不变的描述（Proportional to the principle curvature of the image surface and form a rotationally invariant description of M）。然后，由于M是通过水平和垂直方向的梯度来近似的，他们不是真正的旋转不变。

同样，与Moravec算子一样，我们再来看下面的四张张图：

图中A表示在一个物体的内部或背景上，窗口内的灰度值相对不变，因此该窗口表面上几乎没有曲率，因此M的特征值相对很小；B窗口在一个边缘处，垂直于边缘的地方将有明显很大的曲率，而平行于边缘的地方几乎没什么曲率，因此该形式下M的特征值一个会比较大，另一个较小；C和D对应于角度和离散点，在两个方向都会有很大的曲率，因此，M的特征值都将会很大。假设r1和r2是M的两个特征值，通过上面的分析，可以将一个平面表示为以下三个可区分的区域：

Harris&Stephen提出下面的角点性测度（cornerness measure）

$$C(x,y)=det(M)-k(trace(M))^2\\ det(M)=\lambda_{1}\lambda_{2}=AB-C^2\\ trace(M)=\lambda_{1}+\lambda_{2}=A+B$$

k一般取值04~0.6。

计算步骤

1. 对每一个像素计算自相关矩阵M $$M=\left[\begin{matrix}A&C\\C&B\end{matrix}\right]\\ A=\left( \frac{\partial I}{\partial x} \right)^2\otimes\omega_{i}，B=\left( \frac{\partial I}{\partial x} \right)^2\otimes\omega_{i}，C=\left(\frac{\partial I}{\partial x}\frac{\partial I}{\partial y}\right)\otimes\omega_{i}$$
2. 构造角点性映射图（Construct cornerness map） $$C(x,y)=det(M)-k(trace(M))^2\\ det(M)=\lambda_{1}\lambda_{2}=AB-C^2\\ trace(M)=\lambda_{1}+\lambda_{2}=A+B$$
3. 阈值化，对得到的C(x,y)进行阈值
4. 非极大值抑制

总结

Harris算子针对Moravec算子的不足进行了改进，提高了特征点的检测率以及Repeatability。但是，Harris算子计算量大，对尺度很敏感，不具有尺度不变形；另外Harris对特征点的定位也不是很精确，而且Harris也是各向异性的，对噪声敏感。
两篇参考：
https://blog.csdn.net/kezunhai/article/details/11265167
https://blog.csdn.net/xiaxiazls/article/details/8184124

---

SUSAN算子

SUSAN算子（Smallest Univalue Segment Assimilating Nucleus）是一种高效的边缘和角点检测算子，并且具有结构保留的降噪功能（structure preserving noise reduction )。

SUSAN算子通过用一个圆形模板在图像上移动，一般这个圆形模板的半径是（3.4pixels）的包含37个像素。模板内的每一个像素与中心像素进行比较。

为了介绍和分析的需要，我们首先来看下面这个图：

该图是在一个白色的背景上，有一个深度颜色的区域（dark area），用一个圆形模板在图像上移动，若模板内的像素灰度与模板中心的像素（被称为核Nucleus）灰度值小于一定的阈值，则认为该点与核Nucleus具有相同的灰度，满足该条件的像素组成的区域就称为USAN（Univalue Segment Assimilating Nucleus）。

接下来，我们来分析下上图中的五个圆形模的USAN值。对于上图中的e圆形模板，它完全处于白色的背景中，根据前面对USAN的定义，该模板处的USAN值是最大的；随着模板a到d的移动，USAN值逐渐减少；当圆形模板移动到b处时，其中心位于边缘直线上，此时其USAN值逐渐减少为最大值的一半；而圆形模板运行到角点处a时，此时的USAN值最小。因此通过上面的描述：我们可以推导出：边缘处的点的USAN值小于或等于最大值一半。由此，我们可以得出SUSAN提取边缘和角点算法的基本原理：在边缘或角点处的USAN值最小，可以根据USAN区域的大小来检测边缘、角点等特征的位置和方向信息。

上面都是口头阐述，文字的力量是单薄的，下面我们进入公式阶段。SUSAN算子通过用一个圆形模板在图像上移动，一般这个圆形模板的半径是（3.4pixels）的包含37个像素。模板内的每一个像素与中心像素进行比较，比较方式如下所示：

$$C(\overrightarrow{r},\overrightarrow{r_{0}})= \begin{equation} \left\{ \begin{array}{lr} 1&if\left| I(\overrightarrow{r})-I(\overrightarrow{r_{0}}) \right|\leq t \\ 0&if\left| I(\overrightarrow{r})-I(\overrightarrow{r_{0}}) \right|< t \\ \end{array} \right. \end{equation}$$

其中r0是中心像素，r是掩膜内的其他像素，t是一个像素差异阈值（通常对于对比度比较低的区域，选取较小的t；反之，则t的阈值可以选择大些）。 接着，对上式进行统计，统计方式如下式：

$$n(\overrightarrow{r_{0}})=\sum_{\overrightarrow{r_{0}}}^{}{C(\overrightarrow{r},\overrightarrow{r_{0}})}$$

得到的n值就是USAN的大小。
得到USAN值后，通过阈值化就可以得到初步的边缘响应，公式表示如下：

$$C(\overrightarrow{r},\overrightarrow{r_{0}})= \begin{equation} \left\{ \begin{array}{lr} g-n(\overrightarrow{r_{0}})&if &n(\overrightarrow{r_{0}}) < g \\ 0 && otherwise \end{array} \right. \end{equation}$$

其中，g=3Nmax/4，也即g的取值为USAN最大值的3/4。USAN值越小，边缘的响应就越强。
对得出的边缘响应进行非极大值抑制，就可以得到图像的边缘信息了。

以上完成了SUSAN检测边缘的功能，或许你已经想到了怎么用SUSAN算子来检测角点了。通过上面对a、b、c、d、e等几个圆形模板的USAN值的分析，当模板的中心位于角点处时，USAN的值最小。下面简单叙述下利用SUSAN算子检测角点的步骤：

1. 利用圆形模板遍历图像，计算每点处的USAN值
2. 设置一阈值g，一般取值为1/2(Max(n)， 也即取值为USAN最大值的一半，进行阈值化，得到角点响应
3. 使用非极大值抑制来寻找角点。

通过上面的方式得到的角点，存在很大伪角点。为了去除伪角点，SUSAN算子可以由以下方法实现：①计算USAN区域的重心，然后计算重心和模板中心的距离，如果距离较小则不是正确的角点；②判断USAN区域的重心和模板中心的连线所经过的像素都是否属于USAN区域的像素，如果属于那么这个模板中心的点就是角点。

总结

SUSAN算子是一个原理简单、易于了解的算子。由于其指数基于对周边象素的 灰度比较，完全不涉及梯度的运算，因此其抗噪声能力很强，运算量也比较小；同时，SUSAN算子还是一个各向同性的算子；最后，通过控制参数t和g，可以根据具体情况很容易地对不同对比度、不同形状的图像通过设置恰当的t和g进行控制。比如图像的对比度较大，则可选取较大的t值，而图像的对比度较小，则可选取较小的t值。总之，SUSAN算子是一个非常难得的算子，不仅具有很好的边缘检测性能；而且对角点检测也具有很好的效果。
参考：
http://blog.csdn.net/kezunhai/article/details/11269793
http://users.fmrib.ox.ac.uk/~steve/susan/
http://blog.csdn.net/augusdi/article/details/9012555

---

FAST算子

FAST算子是Rosten等人在SUSAN角点特征检测方法的基础上利用机器学习方法提出的。全名是：Features from Accelerated Segment Test。它计算速度快，可以应用与实时场景中。在FAST特征提出之后，实时计算机视觉应用中特征提取性能才有显著改善。

目前以其高计算效率(computational performance)、高可重复性(high repeatability)成为计算机视觉领域最流行的角点检测方法。
FAST算法包含3个主要步骤：

1. 对固定半径圆上的像素进行分割测试，通过逻辑测试可以去处大量的非特征候选点；
2. 基于分类的角点特征检测，利用ID3 分类器根据16个特征判决候选点是否为角点特征，每个特征的状态为-1，0，1。
3. 利用非极大值抑制进行角点特征的验证。

下面我们对这3个步骤进行分析。

固定半径圆上的像素进行分割测试（Segment Test）

在分割测试中，可以去除大量的非候选角点，这样就可以把可能的角点筛选出来。分割测试是通过对一固定半径的圆形模板的比较和计算进行的，在FAST角点检测算子中，一般是通过半径为3.4 pixel、外围16个像素的圆的作为模板，通过下图我们来具体分析下分割测试。

在上图中，是12点的分割测试角点检测算法的示意图（还有9点的分割测试角点检测，具体参考后面的参考资料）。12点分割测试角点检测算法是在一个图像块上进行，如上图中的左边方块，其中p是中心像素点，12点取的是图上用弧线连接的12个点的像素值都大于或都小于中心像素值，则认为该点处是候选角点（为什么选择12点，因为通过测试，12点的角点检测性能最稳定、速度更快、效果也很好，当然有些文献指出9点的方式也很好）。分割测试是怎么进行的呢？用下面的公式来说，便可一目了然。

$$S_{p\rightarrow x}= \begin{equation} \left\{ \begin{array}{lr} d, I_{p\rightarrow x}<I_p-t & darker\\ s, I_p-t<I_{p\rightarrow x}<I_p+t& similar\\ b,I_p+t<I_{p\rightarrow x} & brighter \end{array} \right. \end{equation}$$

上式中，t是一个阈值（默认取值为10，不同场景取值有差异），Ip表示的是中心像素的像素值，Ip->x表示的是圆形模板中的像素值上式的意思是：当中心像素的像素值Ip小于x处的像素值Ip->x+t时，则该像素属于darker，Sp->x=d，其他两种情况分别表示亮些和相似。这样一个块（圆形）区域就可以分成d、s和b三种类型。这时候只要统计圆形区域中d或b的次数，只要d或b出现的次数大于n（（当是12点分割测试角点检测时，n=12；当是9点时，则n=9），那么该点就被认为是候选角点。

在分割测试步骤中，为了加快速度，其实不需要对这些像素进行逐一的比较。简单来说：首先比较1、5、9、13处点的像素值（也即水平方向和垂直方向上的4个点）与中心像素值的大小，如果这四个点中的像素值有3个或3个以上大于Ip+t或小于Ip-t，那么则认为该处是一个候选角点，否则就不可能是角点。

ID3决策树算法来训练角点检测

通过上式中比较，可以将模板内的像素分成三部分d、s、b，分别记为：Pd，Ps，Pb。因此对于每个Sp->x都属于Pd，Ps，Pb中的一个。另外，令Kp为true，如果p为角点，否则为false。通过ID3算法来选择具有最大信息增益的像素来判断一个像素是否为角点。Kp的熵用下式来计算：

$$c=|\left\{ p|K_{p}=true \right\}|，角点\\ \bar{c}=|\left\{ p|K_{p}=false \right\}|，非角点$$

某一像素的信息增益通过下式来表示：

$$H(P)-H(P_{d})-H(Ps)-H(P_{b})$$

对上述像素依次进行如上处理，选择像素增益最大的像素作为判断角点的依据，生成决策树，从而实现角点的正确分类。

非极大值抑制

在上面的分割测试中，没有计算角点响应函数（Corner Response Function），非极大值抑制无法直接应用于提取的特征。因此，定义一个角点响应函数V，考虑到分割测试的特征以及计算速度的需要，角点响应函数的定义如下：

$$V=\max(\sum_{x\in S_{bright}}^{}{|I_{p\rightarrow x}-I_{p}|-t}，\sum_{x\in S_{dark}}^{}{|I_{p}-I_{p\rightarrow x}|-t})$$

定义了角点响应函数后，就可以采用常规的非极大值抑制方法对非角点进行排除了。

总结

FAST角点检测算法是一种具有高计算效率(computational performance)、高可重复性(high repeatability)特征提取算子，在立体图像匹配、图像配准、目标识别、目标跟踪、场景重构等领域得到了广泛的应用，成为计算机视觉领域最流行的角点检测方法。但是，噪声对该算子的影响比较大，而且阈值t对算子的影响比较大。
参考：
http://blog.csdn.net/kezunhai/article/details/11290749
FAST Corner Detection — Edward Rosten
A Brief History of FAST corner detector
AGAST Corner Detector: faster than FAST and even FAST-ER

---

DoG算子

DoG算子是由Lowe D.G.提出的，对噪声、尺度、仿射变化和旋转等具有很强的鲁棒性，能够提供更丰富的局部特征信息。

作为一个增强算法，DOG可以被用来增加边缘和其他细节的可见性，大部分的边缘算子使用增强高频信号的方法，但是因为随机噪声也是高频信号，很多锐化算子也增强了噪声。DOG算法去除的高频信号中通常包含了随机噪声，所以这种方法是最适合处理那些有高频噪声的图像。这个算法的一个主要缺点就是在调整图像对比度的过程中信息量会减少。

当它被用于图像增强时，DOG算法中两个高斯核的半径之比通常为4:1或5:1。当k设为1.6时，即为高斯拉普拉斯算子的近似。高斯拉普拉斯算子在多尺度多分辨率像片中，用于近似高斯拉普拉斯算子两个高斯核的确切大小决定了两个高斯模糊后的影像间的尺度。

DOG也被用于尺度不变特征变换中的斑点检测。事实上，DOG算法作为两个多元正态分布的差通常总额为零，把它和一个恒定信号进行卷积没有意义。当K约等于1.6时它很好的近似了高斯拉普拉斯变换，当K约等于5时又很好的近似了视网膜上神经节细胞的视野。它可以很好的作为一个实时斑点检测算子和尺度选择算子应用于递归程序。

推导

我们这次从数学推导开始，以LoG作为铺垫，再介绍DoG。

Difference of Gaussian(DOG)是高斯函数的差分。它是可以通过将图像与高斯函数进行卷积得到一幅图像的低通滤波结果，即去噪过程，这里的Gaussian和高斯低通滤波器的高斯一样，是一个函数，即为正态分布函数。同时，它是对Laplacian of Gaussian(LoG)的近似，在某一尺度上的特征检测可以通过对两个相邻高斯尺度空间的图像相减，得到DoG的响应值图像。

Laplace算子对通过图像进行操作实现边缘检测时，对离散点和噪声比较敏感。于是，首先对图像进行高斯暖卷积滤波进行降噪处理，再采用Laplace算子进行边缘检测，就可以提高算子对噪声和离散点的Robust, 这一个过程中Laplacian of Gaussian(LOG)算子就诞生了。

高斯卷积函数定义为：

$$G_{\sigma_{1}}(x,y)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left( -\frac{x^2+y^2}{2\sigma^2} \right)$$

而原始图像与高斯卷积定义为：

$$\Delta[G_{\sigma}(x,y)\ast f(x,y)]=[\Delta G_{\sigma}(x,y)]\ast f(x,y)=LoG\ast f(x,y)$$

因为：

$$\frac{d}{dt}[h(t)\ast f(t)] = \frac{d}{dt}\int_{}^{}f(\tau)h(t-\tau)d\tau=\int_{}^{}f(\tau)\frac{d}{dt}h(t-\tau)d\tau=f(t)\ast \frac{d}{dt}h(t)$$

所以Laplacian of Gaussian(LOG)

$$\Delta G_{\sigma}(x,y)$$

可以通过先对高斯函数进行偏导操作，然后进行卷积求解。公式表示为：

$$\frac{\partial}{\partial x}G_{\sigma}(x,y)=\frac{\partial}{\partial x}e^{-(x^2+y^2)/2\sigma^2}=-\frac{1}{\sigma^2}e^{-(x^2+y^2)/2\sigma^2}$$ $$\frac{\sigma^2}{\sigma^2x}G_{\sigma}(x,y)=\frac{x^2}{\sigma^4}e^{-(x^2+y^2)/2\sigma^2}-\frac{1}{\sigma^2}e^{-(x^2+y^2)/2\sigma^2}=\frac{x^2-\sigma^2}{\sigma^4}e^{-(x^2+y^2)/2\sigma^2}$$

因此，我们可以把LOG核函数定义为：

$$LoG\stackrel{\Delta}{=}\Delta G_{\sigma}(x,y)=\frac{\sigma^2}{\sigma x^2}G_{\sigma}(x,y)+\frac{\sigma^2}{\sigma y^2}G_{\sigma}(x,y)=\frac{x^2+y^2-2\sigma^2}{\sigma^4}e^{-(x^2+y^2)/2\sigma^2}$$ $$\begin{equation} \begin{split} \nabla^2g(x,y)&=\frac{\partial^2g(x,y,\sigma)}{\partial x^2}U_{x}+\frac{\partial^2g(x,y,\sigma)}{\partial y^2}U_{y} \\&= \frac{\partial\nabla g(x,y,\sigma)}{\partial x}U_{x}+\frac{\partial\nabla g(x,y,\sigma)}{\partial y}U_{y} \\&= \left( \frac{x^2}{\sigma^2}-1 \right)\frac{e^{-(x^2+y^2)/2\sigma^2}}{\sigma^2}+\left( \frac{y^2}{\sigma^2}-1 \right)\frac{e^{-(x^2+y^2)/2\sigma^2}}{\sigma^2} \\&= \frac{1}{\sigma^2}\left( \frac{x^2+y^2}{\sigma^2}-2 \right)e^{-(x^2+y^2)/2\sigma^2} \end{split} \end{equation}$$

高斯函数和一级、二阶导数如下图所示：

Laplacian of Gaussian计算可以利用高斯差分来近似，其中差分是由两个高斯滤波与不同变量的卷积结果求得的

$$\sigma\nabla^2g(x,y,\sigma)=\frac{\partial g}{\partial\sigma}\approx\frac{g(x,y,\sigma)-g(x,y,\sigma)}{k\sigma-\sigma}$$

Difference of Gaussian(DOG)是高斯函数的差分。它是可以通过将图像与高斯函数进行卷积得到一幅图像的低通滤波结果，即去噪过程，这里的Gaussian和高斯低通滤波器的高斯一样，是一个函数，即为正态分布函数。同时，它对高斯拉普拉斯LoG(博文LOG算子介绍了实现原理)的近似，在某一尺度上的特征检测可以通过对两个相邻高斯尺度空间的图像相减，得到DoG的响应值图像。

两幅图像的高斯滤波表示为：

$$g_{1}(x,y)=G_{\sigma_{1}}(x,y)*f(x,y) \\ g_{2}(x,y)=G_{\sigma_{2}}(x,y)*f(x,y)$$

最后，将上面滤波得到的两幅图像g1和g2相减

$$g_{1}(x,y)-g_{2}(x,y)=G_{\sigma_{1}}*f(x,y) - G_{\sigma_{2}}*f(x,y)=(G_{\sigma_{1}}-G_{\sigma_{2}})*f(x,y)=DoG*f(x,y)$$

即：可以DOG表示为：

$$DoG\stackrel{\Delta}{=}G_{\sigma_{1}}- G_{\sigma_{2}}=\frac{1}{\sqrt{2\pi}}\left( \frac{1}{\sigma_{1}}e^{-(x^2+y^2)/2\sigma_{1}^2} - \frac{1}{\sigma_{2}}e^{-(x^2+y^2)/2\sigma_{2}^2} \right)$$

在具体图像处理中，就是将两幅图像在不同参数下的高斯滤波结果相减，得到DoG图。

标记红色当前像素点，黄色的圈标记邻接像素点，用这个方式，最多检测相邻尺度的26个像素点。如果它是所有邻接像素点的最大值或最小值点，则标记红色被标记为特征点，如此依次进行，则可以完成图像的特征点提取。

参考

Laplacian of Gaussian
Difference of Gaussian(DOG)
http://www.thinkface.cn/thread-626-1-1.html
http://blog.csdn.net/xiaowei_cqu/article/details/8067881
http://blog.csdn.net/xiaowei_cqu
http://blog.sina.com.cn/s/blog_4bdb170b0101ovif.html
http://blog.csdn.net/abcjennifer/article/details/7639488
http://www.tandfonline.com/doi/abs/10.1080/757582976
http://www.voidcn.com/blog/hh555800/article/p-3597927.html
http://fourier.eng.hmc.edu/e161/lectures/gradient/node9.html