【计算机视觉】2. 特征点检测：Harris, SIFT, SURF, ORB

特征点检测：Harris, SIFT, SURF & ORB

Harris角点检测

Def. [角点（corner point）]
在邻域内的各个方向上灰度变化值足够高的点，是图像边缘曲线上曲率极大值的点。

• [基于灰度图像的角点检测] 包括基于梯度的方法（通过计算边缘的曲率判断角点），基于模板的方法（考虑像素邻域点的灰度变化, 将与邻点亮度对比足够大的点定义为角点），基于模板梯度组合的方法
• [基于二值图像的角点检测] 将二值图像作为单独的检测目标，可使用各种基于灰度图像的角点检测方法
• [基于轮廓曲线的角点检测] 通过角点强度或曲线曲率提取角点

角点检测的基本思想是：使用角点检测算子，对图像的每个像素计算角点响应函数(Corner Response Function )，阈值化角点响应函数，根据实际情况选择阈值，对阈值化的角点响应函数进行非极大值抑制，并获取非零点作为角点。通过一个小的滑动窗口在邻域检测角点
在任意方向上移动窗口，若窗口内的灰度值都有剧烈的变化，则窗口的中心就是角点。定义角点响应函数 E(u,v)=\sum_{(x,y)} w(x,y) [I(x+u,y+v)-I(x,y)]^2

其中w为窗函数（window function），I为图像梯度，那么E就表示了灰度变化的剧烈程度。1977年，Moravec最先提出了如下的角点检测方法：

1. 对于原始图像，取偏移量(Δx,Δy)为(1,0),(1,1),(0,1),(-1,1)，分别计算每一像素点(xi,yi)的灰度变化
2. 对于每一像素点(xi,yi)，计算角点响应函数R(xi,yi)=min E
3. 设定阈值T，将角点响应函数R(xi,yi)中低于T的值设为0
4. 在窗口范围内进行非极大值抑制：遍历角点响应函数，若某个像素的角点响应函数在窗口内不是最大，该像素置0
5. 选择非零点作为角点检测结果

Moravec角点检测的缺点

• 二值的窗口函数导致角点响应函数不够光滑
• 只在四个方向上计算灰度值变化，导致角点响应函数在多处都有较大响应
• 对于每个点只考虑E的最小值，导致算法对边缘有很强的反应

1988年，Harris和Plessey对Moravec的方法进行了改进，提出了经典的Harris角点检测算法。Harris首先将Moravec算法中的窗口函数由阶跃函数改为二维高斯函数，并通过泰勒展开考察微小移动，也就是说，如果要求E的最大值以明确角点，就可以令$u,v \rightarrow 0$，对E做泰勒展开，得

E(u,v)=(u,v)M\begin{pmatrix} u\\ v \end{pmatrix}\\ \text{where }M=\sum_{(x,y)}w(x,y)\begin{pmatrix} I_X^2 & I_X I_Y \\ I_X I_Y & I_Y^2 \end{pmatrix}=\begin{pmatrix} \sum_W I_X^2 & \sum_W I_X I_Y \\ \sum_W I_X I_Y & \sum_W I_Y^2 \end{pmatrix}

记 M=\begin{pmatrix}A & B\\ B & C\end{pmatrix} ，则上式可以写成 Au^2+2Buv+Cv^2=E 的形式，这表示了一个椭圆，自相关矩阵M描述了图像局部区域的灰度变化趋势，可以通过椭圆的形状来判定角点。

Proposition. 对于椭圆 Ax^2+2Bxy+Cy^2=1 ，设其半长轴和半短轴分别为a,b，那么 \frac{1}{\sqrt{a}},\frac{1}{\sqrt{b}} 是矩阵M的特征值。

Proof. 椭圆上任意一点到原点（也就是椭圆的中心）的距离平方 d^2=x^2+y^2 ，条件为 Ax^2+2Bxy+Cy^2-1=0 ，那么拉格朗日函数 L=x^2+y^2+\lambda(Ax^2+2Bxy+Cy^2-1) 于是
\begin{cases} F_x'=(2-2A\lambda)x-2B\lambda y=0\\ F_y'=-2B\lambda x+(2-2C\lambda)y=0 \end{cases}
则
[(1-A\lambda)x^2-B\lambda xy]+[-B\lambda xy+(1-C\lambda)y^2]=0
又可以改写为
\begin{pmatrix} x & y \end{pmatrix}\begin{pmatrix} 1-A\lambda & -B\lambda \\ -B\lambda & 1-C\lambda \end{pmatrix}\begin{pmatrix} x\\y \end{pmatrix}=0
由于（x,y）不是零向量，所以 \begin{vmatrix}1-A\lambda & -B\lambda \\ -B\lambda & 1-C\lambda\end{vmatrix}=0 ，即1/λ是矩阵 M 的特征值，所以
\begin{cases} \frac{1}{\lambda_1}+\frac{1}{\lambda_2}=A+C=\text{tr} M\\ \frac{1}{\lambda_1 \lambda_2}=AC-B^2=\det M \end{cases}
又

d^2=x^2+y^2=\lambda(Ax^2+2Bxy+Cy^2)=\lambda\Rightarrow d=\sqrt{\lambda}

这样λ1，λ2就分别对应d^2的两个最值，也就是a^2和b^2，所以， \frac{1}{\sqrt{a}},\frac{1}{\sqrt{b}} 是矩阵M的特征值。由此我们还可以得到，椭圆面积 S=\pi ab=\frac{\pi}{\sqrt{AC-B^2}}
q.e.d.

得到M的特征值有什么用呢？若用奇异值分解的观点看这个问题。由于Jacobian矩阵 J=\begin{pmatrix} I_X \\ I_Y\end{pmatrix} ，故M=JJ^T，这样M的特征值开根号后就是J的奇异值，因此M的特征值就可以体现IX和IY的相对大小。

现在我们需要定义角点响应函数，进一步进行区分，令（一般k=0.04 ~ 0.06）
R=\text{det}(M)-k \text{tr}^2 (M)
其中
\begin{cases} \text{det}⁡(M)=\lambda_1 \lambda_2=\sum_W I_X^2 \cdot \sum_W I_Y^2 - (\sum_W I_X I_Y)^2\\ \text{tr}(M)=\lambda_1+\lambda_2 = \sum_W I_X^2+I_Y^2 \end{cases}
则定性判断方法为

• [边缘] \lambda_2\geq \lambda_1 或\lambda_2\leq \lambda_1
• [角点] \lambda_1,\lambda_2 都很大，E在各个方向显著变化
• [平滑区域] \lambda_1,\lambda_2 都很小，$E$在各个方向基本不变

定量判断方法为

• [边缘] R<<0
• [角点] R>>0
• [平滑区域] |R|很小

一般增大k的值，将减小角点响应值R，降低角点检测的灵性，减少被检测角点的数量；减小k值，将增大角点响应值R，增加角点检测的灵敏性，增加被检测角点的数量。

下面我们介绍如何利用Harris角点特征进行特征点匹配。

对两幅图像进行特征匹配的过程是：

1. 建立图像的特征点数据库每个特征点的数据结构，包括：位置坐标、尺度、方向、特征向量，
2. 为新图像的每个特征点在数据库中逐个匹配，根据特征向量的欧氏距离在数据库中寻找其最近邻和次近邻特征点，若(最近邻距离/次近邻距离)大于某一阙值，则特征匹配成功。

Harris 角点检测器仅仅能够检测出图像中的兴趣点，但是没有给出通过比较图像间的兴趣点来寻找匹配角点的方法。我们需要在每个点上加入描述子信息，并给出一个比较这些描述子的方法。

兴趣点描述子是分配给兴趣点的一个向量，描述该点附近的图像的表观信息。描述子越好，寻找到的对应点越好。我们用对应点或者点的对应来描述相同物体和场景点在不同图像上形成的像素点。Harris 角点的描述子通常是由周围图像像素块的灰度值，以及用于比较的归一化互相关矩阵构成的。图像的像素块由以该像素点为中心的周围矩形部分图像构成。通常，两个（相同大小）像素块 I1(x) 和 I2(x) 的归一化互相关矩阵定义为：
ncc(I_1,I_2 )=\frac{1}{n-1} \sum_x \frac{[I_1 (x)-\mu_1]}{\sigma_1} \frac{[I_2 (x)-\mu_2]}{\sigma_2}

Harris角点的优点

• 计算简单
• 提取的点特征均匀且合理
• 稳定：Harris算子对图像旋转、亮度变化、噪声影响和视点变换不敏感

Harris 算子的局限性

• 对尺度很敏感，不具有尺度不变性
• 提取的角点精度是像素级的
• 需要设计角点匹配算法

1994年，Shi和Tomasi对Harris角点检测进行了改进，即定义角点响应函数R=min⁡ (λ1，λ2)，该改进是计算适合跟踪的优质特征（good features to track），使得特征分布更均匀，其检测精度是亚像素级别的。

SIFT算法预备知识：尺度空间与图像金字塔

1999年Lowe提出了SIFT(Scale-invariant feature transform，尺度不变性特征变换)特征检测算法，并于2003年对其完善总结；2006年，Bay等人对SIFT算法进行了改进，提升其效率，提出了SURF(Speeded Up Robust Features，加速鲁棒性特征)算法。

SIFT特征检测算法的特点：

• SIFT特征是图像的局部特征，对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性
• 信息量丰富，适用于在海量特征数据库中进行匹配
• 多量性，少数物体也可以产生大量SIFT特征
• 高速性，经优化的SIFT匹配算法甚至可以达到实时性

SIFT特征检测的步骤：

• 检测尺度空间的极值点
• 精确定位特征点(Keypoint)
• 设定特征点的方向参数
• 生成特征点的描述子（128维向量）

如果一张照片的像素数目被不断压缩，或者观察者距离照片越来越远，它将逐渐变得模糊，而这个导致照片的呈现内容发生变化的连续的自变量就可以称为尺度。对物体观察的尺度不同，物体呈现的方式也不同。对计算机视觉而言，无法预知某种尺度的物体结构是否有意义，因此有必要将所有尺度的结构表示出来。以下关于尺度空间理论的部分，也可以参考：http://www.cnblogs.com/ronny/p/3886013.html

利用图像金字塔的方法，我们可以得到一系列大小不一的图像，由大到小，从下到上构成的塔状模型。假设高斯金字塔的第l 层图像为Gl，则有：
G_l(i,j)=\sum^{2}_{m=-2}\sum_{n=-2}^2\omega(m,n)G_{l-1}(2i+m,2j+n)\\ (1\le l \le N,0 \le i \le R_l,0 \le j \le C_l)
式中，N为高斯金字塔层数；Rl和Gl分别为高斯金字塔第l层的行数和列数；ω(m,n)是一个二维可分离的5× 5窗口函数，表达式为：
\omega=\frac{1}{256} \begin{pmatrix} 1 & 4 & 6 & 4 & 1\\ 4 & 16 & 24 & 16 & 4\\ 6 & 24 &36 & 24 & 6\\ 4 & 16 & 24 & 16 & 4\\ 1 & 4 & 6 & 4 & 1 \\ \end{pmatrix} =\frac{1}{16}\begin{pmatrix} 1 & 4 & 6 & 4 &1 \end{pmatrix}\times\frac{1}{16}\begin{pmatrix} 1 \\ 4 \\ 6 \\ 4 \\1 \end{pmatrix}
写成上面的形式是为了说明，二维窗口的卷积算子，可以写成两个方向上的1维卷积核（二项核）的乘积。上面卷积形式的公式实际上完成了2个步骤：1）高斯模糊；2）降维。按上述步骤生成的G0,G1,...,GN就构成了图像的高斯金字塔，其中$G0$为金字塔的底层(与原图像相同)，$GN$为金字塔的顶层。可见高斯金字塔的当前层图像是对其前一层图像先进行高斯低通滤波，然后做隔行和隔列的降采样(去除偶数行与偶数列)而生成的。

将Gl进行内插（这里内插用的不是双线性而是用的与降维时相同的滤波核）得到放大图像 G_l^\ast ，使 G_l^\ast 的尺寸与 G_{l-1} 的尺寸相同，表示为：
G_l^{\ast}(i,j)=4\sum_{m=-2}^2\sum_{n=-2}^2\omega(m,n)G_l(\frac{i+m}{2},\frac{j+n}{2})
其中 0 \le l \le N,0 \le i \le R_l,0 \le j \le G_l ，上面的系数4，是因为每次能参与加权的项，权值和为4/256，这个与我们用的ω的值有关。式中，
G_l(\frac{i+m}{2},\frac{j+n}{2})=\begin{cases} G_l(\frac{i+m}{2},\frac{j+n}{2}), & \frac{i+m}{2},\frac{j+n}{2}\in \mathbb{Z} \\ 0, & \text{others} \end{cases}
令
\begin{cases} LP_l=G_l-G_{l+1}^{\ast}, & 0 \le l \le N \\ LP_N=G_N, & l=N \end{cases}
式中，N为拉普拉斯金字塔顶层号，LPl是拉普拉斯金字塔分解的第$l$层图像。由LP0,LP1,\cdots,LPl,...,LPN构成的金字塔即为拉普拉斯金字塔。它的每一层图像是高斯金字塔本层图像与其高一级的图像经内插放大后图像的差，此过程相当于带通滤波，因此拉普拉斯金字塔又称为带通金字塔分解。

图像的金字塔化能高效地（计算效率也较高）对图像进行多尺度的表达，但它缺乏坚实的理论基础，不能分析图像中物体的各种尺度。

我们知道，信号的尺度空间刚提出是就是通过一系列单参数、宽度递增的高斯滤波器将原始信号滤波得到到组低频信号。那么一个很明显的疑问是：除了高斯滤波之外，其他带有参数t的低通滤波器是否也可以用来生成一个尺度空间。后来Koenerink、Lindeberg、Florack等人用精确的数学形式通过不同的途径都证明了高斯核是实现尺度变换的唯一变换核。

虽然很多研究者从可分性、旋转不变性、因果性等特性推出高斯滤波器是建立线性尺度空间的最优滤波器。然后在数字图像处理中，需要对核函数进行采样，离散的高斯函数并不满足连续高斯函数的的一些优良的性质。所以后来出现了一些非线性的滤波器组来建立尺度空间，如B样条核函数。

使用高斯滤波器对图像进行尺度空间金塔塔图的构建，让这个尺度空间具有下面的性质：

• 加权平均和有限孔径效应
• 信号在尺度t上的表达可以看成是原信号在空间上的一系列加权平均，权重就是具有不同尺度参数的高斯核。
• 信号在尺度t上的表达也对应于用一个无方向性的孔径函数（特征长度为 \sigma=\sqrt{t} ）来观测信号的结果。这时候信号中特征长度小于σ的精细结构会被抑制，理解为一维信号上小于σ的波动会被平滑掉。
• 层叠平滑，也叫高斯核族的半群（Semi-Group）性质：两个高斯核的卷积等同于另外一个不同核参数的高斯核卷积。 g(\mu,\sigma_1)\ast g(\mu,\sigma_2)=g(\mu,\sqrt{\sigma_1^2+\sigma_2^2}) 这个性质的意思就是说不同的高斯核对图像的平滑是连续的。

• 局部极值递性
• 这个特征可以从人眼的视觉原理去理解，人在看一件物体时，离得越远，物体的细节看到的越少，细节特征是在减少的。对生理学的研究中发现，哺乳动物的视网膜和视觉皮层的感受区域可以很好地用4阶以内的高斯微分来建模。
• 高斯核对图像进行滤波具有压制局部细节的性质。
• 尺度伸缩不变性，对原来的信号加一个变换函数，对变换后的信号再进行高斯核的尺度空间生成，新的信号的极值点等特征是不变的。

SIFT算法

令原图像为金字塔的第一层，每次降采样所得到的新图像为金字塔的一层(每层一张图像)，每个金字塔共n层。金字塔的层数根据图像的原始大小和塔顶图像的大小共同决定，其计算公式如下：
n=\log_2⁡ [\min (M,N)] –t
其中M，N为原图像的长和宽，t为塔顶图像的最小维数的对数值。接下来，我们通过尺度空间理论模拟图像数据的多尺度特征 ，以高斯函数为实现尺度变换的唯一线性核，则二维图像I(x, y)的尺度空间
L(x,y,\sigma)=I(x,y)*g(x,y,\sigma)\\ \text{where } g(x,y,\sigma)=\frac{1}{2\pi\sigma^2 } \exp⁡ \left(-\frac{x^2+y^2}{2\sigma^2 }\right)
根据高斯函数的性质，我们知道高斯窗的宽度约为6σ，即在[-3σ，3σ]外的区间外，信号通过高斯滤波器的响应为0，因此，高斯函数不同的标准差就可以作为不同尺度的量度，事实上对高斯函数进行傅里叶变换后有 \mathscr{F}[g(x,k\sigma)]=k\mathscr{F}[g(x,\sigma)]\quad \mathscr{F}[g(x,y,k\sigma)]=k^2\mathscr{F}[g(x,y,\sigma)] 即一维高斯核是线性核，二维高斯核使得 x\leftarrow kx,y\leftarrow ky ，所以仍是线性的。这样就可以定义高斯差分（difference of Gaussian, DoG）算子
D(x,y,\sigma,k)=L(x,y,k\sigma)-L(x,y,\sigma) \approx (k-1) \nabla^2 h
其中 \nabla^2 h 表示高斯拉普拉斯算子（LoG），正则化的高斯-拉普拉斯变换

\nabla^2_{norm}=\sigma^2\nabla^2g=\sigma^2(\frac{\partial^2g}{\partial x^2}+\frac{\partial^2g}{\partial y^2}) = -\frac{1}{2\pi\sigma^2}[1-\frac{x^2+y^2}{\sigma^2}]\cdot \exp(-\frac{x^2+y^2}{2\sigma^2})

求 \frac{\partial(\nabla^2_{norm})}{\partial\sigma} = 0\Rightarrow r^2-2\sigma^2=0

图像与某一个二维函数进行卷积运算实际就是求取图像与这一函数的相似性。同理，图像与高斯拉普拉斯函数的卷积实际就是求取图像与高斯拉普拉斯函数的相似性。当图像中的斑点尺寸与高斯拉普拉斯函数的形状趋近一致时，图像的拉普拉斯响应达到最大。

从概率的角度解释为：假设原图像是一个与位置有关的随机变量X的密度函数，而LOG为随机变量Y的密度函数，则随机变量X+Y的密度分布函数即为两个函数的卷积形式。如果想让X+Y能取到最大值，则X与Y能保持步调一致最好，即X上升时，Y也上升，X最大时，Y也最大。实际上 \nabla^2[G(x,y)*f(x,y)] = \nabla^2[G(x,y)]*f(x,y)

每一次采样返回的结果，就像音乐上的一个八度（octave），为了让尺度体现其连续性，高斯金字塔在简单降采样的基础上加上了高斯滤波。将图像金字塔每层的一张图像使用不同参数做高斯模糊，使得金字塔的每层含有多张高斯模糊图像，将金字塔每层多张图像合称为八度(Octave)，金字塔每层只有一组图像，组数和金字塔层数相等，每组含有多张(也叫层，Interval)图像。另外，降采样时，高斯金字塔上一组图像的初始图像(底层图像)是由前一组图像的倒数第三张图像隔点采样得到的。在下一个八度，高斯算子的标准差扩大一倍。

在高斯差分尺度空间检测局部极大或极小值，检测点与其同尺度的8个相邻点、上下相邻尺度对应的9 × 2个点进行比较，以确保在尺度空间和二维图像空间都检测到极值点，极值点的位置可以通过对高斯差分算子求一阶导数得到。离散空间的极值点并不是真正的极值点，利用已知的离散空间点插值得到的连续空间极值点的方法叫做亚像素插值（Sub-pixel Interpolation）。为了提高关键点的稳定性，需要对尺度空间DoG函数进行曲线拟合。利用DoG函数在尺度空间的Taylor展开式(拟合函数)为：
D(X)=D+\frac{\partial D^T}{\partial X} X+\frac{1}{2} X^T \frac{\partial^2 D}{\partial X^2} X \\ \text{where } X=[x,y,\sigma,k]^T
求导并让方程等于零，可以得到极值点的偏移量为：
\hat X= -\frac{\partial^2 D^{-1}}{\partial X^2 } \frac{\partial D}{\partial X}
对于这样得到的极值点，首先要舍去高斯差分算子绝对值很小（一般指|D(\hat X)|<0.03）的点，通过主曲率分析去除边缘响应过大的极值点，即计算差分图像D的Hessian矩阵：
H= \begin{pmatrix} D_{XX} & D_{XY} \\ D_{XY} & D_{YY} \end{pmatrix}
保留满足如下条件的极值点（一般r=10）
\frac{\text{tr}^2 H}{\text{det} H} < \frac{(r+1)^2}{r}
现在对L(x,y,σ)进行差分，得到邻域梯度方向直方图（这里常用Roberts算子计算），即
\begin{gather} m(x,y)=\sqrt{(L(x+1,y)-L(x-1,y))^2+(L(x,y+1)-L(x,y-1))^2}\\ \theta(x,y)=\arctan \frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)} \end{gather}
进而确定特征的方向参数：

• 梯度方向直方图的范围：[0, 360]；
• 特征点的主方向：直方图的峰值 ；
• 特征点的辅方向：直方图的次峰值，能量超过主峰值能量的80\%；

当图像发生旋转时，特征点的主方向相对于像素的梯度方向不变；将多幅待匹配的图像都旋转到令特征点方向为0的位置再匹配，使特征具有旋转不变性。

现在我们将坐标轴方向旋转为特征点的方向，以特征点为中心取窗口，通过高斯加权增强特征点附近像素梯度方向信息的贡献，即在4 × 4的小块上计算梯度方向直方图（ 取8个方向），计算梯度方向累加值，形成种子点，构成4× 4 × 8= 128维特征向量。SIFT特征对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性。

对描述子的进一步处理

• [去除光照影响] 将特征向量的长度归一化
• [减弱扭曲影响] 仿射变换

关于SIFT算法的几个问题

本小节主要参考：

http://www.cnblogs.com/ronny/p/4028776.html#3962356

第一个问题：第一组第一层图像的生成

这里要分两种情况：其一是把第一组的索引定为0；其二是把第一组的索引定为-1；我们先考虑第一组索引为0的情况，我们知道第一组第一层的图像是由原图像与 \sigma_0
（一般设置为1.6）的高斯滤波器卷积生成，为了图像反走样的需要，通常假设输入图像是经过高斯平滑处理的，其值为 \sigma_n=0.5 ，即半个像元。

也就是说我们采集到的图像I(x,y)，已经被 \sigma=\sigma_n=0.5 的高斯滤波器平滑过了。所以我们不能直接对I(x,y)用 \sigma_0 的高斯滤波器平滑，而应该用 \sigma=\sigma_{20}−\sigma_{2n} 的高斯滤波器去平滑I(x,y)，即
\text{FirstLayer}(x,y) = I_s(x,y)* G(x,y,\sqrt{\sigma_0^2 - (2\sigma_n)^2})
其中FirstLayer(x,y)表示整个尺度空间为第1组第1层的图像。

现在我们来考虑把第一组的索引定为-1的情况。那么首先第一个问题便是为什么要把索引定为-1；如果索引为0，整个尺度空间的第1组的第1层图像已经是由原图像模糊生成的了，那么也就是说已经丢失了细节信息，那么原图像我们完全没有利用上。基于这种考虑，我们先将图像放大2倍，这样原图像的细节就隐藏在了其中。

由上面一种情况分析，我们已经知识了I(x,y)看成是已经被 \sigma_n=0.5 模糊过的图像，那么将I(x,y)放大2倍后得到Is(x,y)，则可以看为是被 2\sigma_n=1 的高斯核模糊过的图像。那么由Is生成第1组第1层的图像用的高斯滤波器的 \sigma = \sqrt{\sigma_0^2 -(2 \sigma_n)^2} ，即
\text{FirstLayer}(x,y) = I_s(x,y)\otimes G(x,y,\sqrt{\sigma_0^2 - (2\sigma_n)^2})

第二个问题：尺度空间生成了多少幅图像

我们知道S是我们最终构建出来的用来寻找特征点的高斯差分图像，而特征点的寻找需要查找的是空间局部极小值，即在某一层上查找局部极值点的时候需要用到上一层与下一层的高斯差分图像，所以如果我们需要查找S层的特征点，需要S+2层高斯差分图像，然后查找其中的第2层到第S+1层。

而每一个高斯差分图像G(x,y,σ )都需要两幅尺度空间的图像L(x,y,kσ )与L(x,y,σ )进行差分生成，这里假设S =3，则我们需要的高斯差分图像有S+2 = 5张，分别为 G(x,y,\sigma ),G(x,y,k\sigma ),G(x,y,k^2\sigma ),G(x,y,k^3\sigma ),G(x,y,k^4\sigma ) 其中的 G(x,y,k\sigma ),G(x,y,k^2\sigma ),G(x,y,k^3\sigma ) 这三张图像是我们用来查找局部极值点的图像。那么我们则需要S+3 = 6张尺度空间图像来生成上面那些高斯差分图像，它们分别为： L(x,y,\sigma ),L(x,y,k\sigma ),L(x,y,k^2\sigma ),L(x,y,k^3\sigma ),L(x,y,k^4\sigma )，L(x,y,k^5\sigma ) 从上面的分析，我们知道对于尺度空间来说，我们一共需要S+3层图像来构建出来S+2层高斯差分图像。所以，如果整个尺度空间一共有n组，每组有S+3层图像，就共有n(S+3)张尺度图像。

第三个问题：为什么取上一张的倒数第3张图像隔行采样后作为下一组的第一张图像？

对于尺度空间第o组，第s层的图像，它的尺度为 \sigma = \sigma_o k^{o+s/S},k =1/2,o\in[0,1,2,\dots,\text{nOctave}-1],s\in[0,1,2,\dots,S+2] 那么我们从第0组开始，看它各层的尺度。

• 第0组 \sigma_0 \to 2^{1/3}\sigma_0\to 2^{2/3}\sigma_0\to 2^{3/3}\sigma_0\to 2^{4/3}\sigma_0\to 2^{5/3}\sigma_0
• 第1组 2\sigma_0\to 2*2^{1/3}\sigma_0\to 2*2^{2/3}\sigma_0\to 2*2^{3/3}\sigma_0\to 2*2^{4/3}\sigma_0\to 2*2^{5/3}\sigma_0
• ……

可以看出，第1组的第0层图像恰好与第0组的倒数第三幅图像一致，尺度都为 2\sigma_0 ，所以我们不需要再根据原图来重新卷积生成每组的第0张图像，只需采用上一层的倒数第3张来降采样即可。

我们也可以继续分析，第0组尺度空间得到的高斯差分图像的尺度为： \sigma_o\to 2^{1/3}\sigma_0\to 2^{2/3}\sigma_0\to 2^{3/3}\sigma_0\to 2^{4/3}\sigma_0 而第1组尺度空间得到的高斯差分图像的尺度为 2\sigma_0\to 2*2^{1/3}\sigma_0\to 2*2^{2/3}\sigma_0\to 2*2^{3/3}\sigma_0\to 2*2^{4/3}\sigma_0 ：如果我们把它们的中间三项取出来拼在一起，则尺度为： 2^{1/3}\sigma_0\to 2^{2/3}\sigma_0\to 2^{3/3}\sigma_0\to 2*2^{1/3}\sigma_0\to 2*2^{2/3}\sigma_0\to 2*2^{3/3}\sigma_0 正好连续。这一效果带来的直接的好处是在尺度空间的极值点确定过程中，我们不会漏掉任何一个尺度上的极值点，而是能够综合考虑量化的尺度因子。

第四个问题：如何用第i-1层的图像生成第i层的图像？

尺度空间里的每一层的图像（除了第1层）都是由其前面一层的图像和一个相对σ的高斯滤波器卷积生成，而不是由原图和对应尺度的高斯滤波器生成的，这一方面是因为我前面提到的不存在所谓意思上的“原图”，我们的输入图像I(x,y)已经是尺度为σ=0.5的图像了。另一方面是由于如果用原图计算，那么相邻两层之间相差的尺度实际上非常小，这样会造成在做高斯差分图像的时候，大部分值都趋近于0，以致于后面我们很难检测到特征点。

基于上面两点原因，所以对于每一组的第i+1层的图像，都是由第i层的图像和一个相对尺度的高斯滤波器卷积生成。

我们首先考虑第0组，它们的第i+1层图像与第i层图像之间的相对尺度为 \text{SigmaDiff}_i=\sqrt{(\sigma_0k^{i+1})^2 – (\sigma_0 k^i)^2} 为了保持尺度的连续性，后面的每一组都用这样一样相对尺度（SIFT实际代码里是这样做的）。这里有一个猜测，比如说尺度为2σ0的这一组，第i层与第i+1层之间的相对尺度计算的结果应该为 \sqrt{(2\sigma_0k^{i+1})^2 – (2\sigma_0 k^i)^2} = 2\cdot \text{SigmaDiff}_i 可是，依然可以用 \text{SigmaDiff}{}_i ，是因为这一层被降维了。

第五个问题：离散空间得到的极值点与连续空间的极值点之间的差别

为了寻找尺度空间的极值点，每一个采样点要和它所有的相邻点进行比较，看其是否比它的图像域和尺度域的相邻点大或小。对于其中的任意一个检测点都要和它同尺度的8个相邻点和上下相邻尺度对应的9× 2个点共26个点比较，以确保在尺度空间和二维图像位置空间都检测到极值点。也就是，比较是在一个3× 3× 3的立方体内进行。

搜索过程从每组的第二层开始，以第二层为当前层，对第二层的DoG图像中的每个点取一个3× 3× 3的立方体，立方体上下层为第一层与第三层。这样，搜索得到的极值点既有位置坐标（DoG的图像坐标），又有空间尺度坐标（层坐标）。当第二层搜索完成后，再以第三层作为当前层，其过程与第二层的搜索类似。当S=3时，每组里面要搜索3层。

极值点的搜索是在离散空间中进行的，检测到的极值点并不是真正意义上的极值点。利用已知的离散空间点插值到连续空间极值点的方法，就是子像元插值方法。

SURF算法

SURF特征(Speeded Up Robust Features，加速鲁棒性特征)是对SIFT特征的进一步优化，基于Hessian矩阵构造金字塔尺度空间，利用箱式滤波器（box filter）简化二维高斯滤波，不需要再进行降采样；通过Haar小波特征设定特征点主方向，这样构建的特征描述子就是64维的；surf构造的金字塔图像与sift有很大不同，就是因为这些不同才加快了其检测的速度。Sift采用的是DOG图像，而surf采用的是Hessian矩阵行列式近似值图像，也写作DOH算子。
H= \begin{pmatrix} L_{XX} & L_{XY} \\ L_{XY} & L_{YY} \end{pmatrix}\Rightarrow \text{det}H=L_{XX} L_{YY}-0.9 L_{XY}^2
由于求Hessian时要先高斯平滑，然后求二阶导数，这在离散的像素点是用模板卷积形成的，这2种操作合在一起用一个模板代替就可以了。相比于sift算法的高斯金字塔构造过程，surf算法速度有所提高。在sift算法中，每一组（octave）的图像大小是不一样的，下一组是上一组图像的降采样（1/4大小）；在每一组里面的几幅图像中，他们的大小是一样的，不同的是他们采用的尺度$\sigma$不同。而且在模糊的过程中，他们的高斯模板大小总是不变的，只是尺度$\sigma$改变。对于surf算法，图像的大小总是不变的，改变的只是高斯模糊模板的尺寸，当然，尺度也是在改变的，但不需要降采样过程，节省时间。

为了保证旋转不变性，在SURF中，不统计其梯度直方图，而是统计特征点领域内的Harr小波特征。即以特征点为中心，计算半径为6s(S为特征点所在的尺度值)的邻域内，统计60度扇形内所有点在x(水平)和y(垂直)方向的Haar小波响应总和(Haar小波边长取4s)，并给这些响应值赋高斯权重系数，使得靠近特征点的响应贡献大，而远离特征点的响应贡献小，然后60度范围内的响应相加以形成新的矢量，遍历整个圆形区域，选择最长矢量的方向为该特征点的主方向。在特征点周围取正方形框，方框的方向为特征点主方向，把方框分为16个子区域，在每个子区域统计水平方向和垂直方向的haar小波特征，在每个子区域计算haar小波特征的水平方向值之和，水平方向绝对值之和、垂直方向值之和、垂直方向绝对值之和，构成16× 4=64维特征向量

在完成采集后，还需要建立图像的特征点数据库，每个特征点的数据结构包括：位置坐标、尺度、方向、特征向量（128或64维）；为新图像的每个特征点在数据库中逐个匹配，即根据特征向量的欧氏距离在数据库中寻找其最近邻和次近邻特征点，若最近邻距离或次近邻距离大于某一阙值，则特征匹配成功。

综上所述，可知SURF采用Henssian矩阵获取图像局部最值还是十分稳定的，但是在求主方向阶段太过于依赖局部区域像素的梯度方向，有可能使得找到的主方向不准确，后面的特征向量提取以及匹配都严重依赖于主方向，即使不大偏差角度也可以造成后面特征匹配的放大误差，从而匹配不成功；另外图像金字塔的层取得不足够紧密也会使得尺度有误差，后面的特征向量提取同样依赖相应的尺度，在这个问题上我们只能采用折中解决方法：取适量的层然后进行插值。

SIFT特征与SURF特征的比较：

• [构建图像金字塔] SIFT特征利用不同尺寸的图像与高斯差分滤波器卷积；SURF特征利用原图片与不同尺寸的方框滤波器卷积。
• [特征描述子] SIFT特征有4×4×8=128维描述子，SURF特征有4×4×4=64维描述子
• [特征点检测方法] SIFT特征先进行非极大抑制，再去除低对比度的点，再通过Hessian矩阵去除边缘响应过大的点；SURF特征先利用Hessian矩阵确定候选点，然后进行非极大抑制
• [特征点主方向] SIFT特征在正方形区域内统计梯度幅值的直方图，直方图最大值对应主方向，可以有多个主方向；SURF特征在圆形区域内计算各个扇形范围内x、y方向的haar小波响应，模最大的扇形方向作为主方向

ORB特征

本节主要参考了：

https://blog.csdn.net/zouzoupaopao229/article/details/52625678

ORB特征是将FAST特征点的检测方法与BRIEF特征描述子结合起来，并在它们原来的基础上做了改进与优化。ORB算法的速度大约是SIFT的100倍，是SURF的10倍。

ORB（Oriented FAST and Rotated BRIEF）是一种快速特征点提取和描述的算法。这个算法是由Ethan Rublee, Vincent Rabaud, Kurt Konolige以及Gary R.Bradski在2011年一篇名为“ORB：An Efficient Alternative to SIFT or SURF”的文章中提出。ORB算法分为两部分，分别是特征点提取和特征点描述。特征提取是由FAST（Features from Accelerated Segment Test）算法发展来的，特征点描述是根据BRIEF（Binary Robust IndependentElementary Features）特征描述算法改进的。

ORB使用ID3算法训练一个决策树，将特征点圆周上的16个像素输入决策树中，以此来筛选出最优的FAST特征点。

接着，非极大值抑制去除局部较密集特征点。使用非极大值抑制算法去除临近位置多个特征点的问题。为每一个特征点计算出其响应大小。计算方式是特征点P和其周围16个特征点偏差的绝对值和。在比较临近的特征点中，保留响应值较大的特征点，删除其余的特征点。

现在我们可以建立金字塔，来实现特征点的多尺度不变性。设置一个比例因子scaleFactor（opencv默认为1.2）和金字塔的层数nlevels（opencv默认为8）。将原图像按比例因子缩小成nlevels幅图像。缩放后的图像为：I’= I/scaleFactork(k=1,2,…, nlevels)。nlevels幅不同比例的图像提取特征点总和作为这幅图像的oFAST特征点。

ORB算法提出使用矩（moment）法来确定FAST特征点的方向。也就是说通过矩来计算特征点以r为半径范围内的质心，特征点坐标到质心形成一个向量作为该特征点的方向。矩定义如下：
m_{pq}=\sum x^p y^q I(x,y)

其中，I(x,y)为图像灰度表达式。该矩的质心为 (m_{10}/m_{00},m_{01}/m_{00})

假设角点坐标为O，则向量的角度即为该特征点的方向。计算公式 \theta=\arctan (m_{10}/m_{01})

rBRIEF特征描述是在BRIEF特征描述的基础上加入旋转因子改进的。BRIEF算法计算出来的是一个二进制串的特征描述符。它是在一个特征点的邻域内，选择n对像素点pi、qi（i=1,2,…,n）。然后比较每个点对的灰度值的大小。如果I(pi)> I(qi)，则生成二进制串中的1，否则为0。所有的点对都进行比较，则生成长度为n的二进制串。一般n取128、256或512，opencv默认为256。另外，值得注意的是为了增加特征描述符的抗噪性，算法首先需要对图像进行高斯平滑处理。

在旋转不是非常厉害的图像里，用BRIEF生成的描述子的匹配质量非常高，作者测试的大多数情况中都超越了SURF。但在旋转大于30°后，BRIEF的匹配率快速降到0左右。BRIEF的耗时非常短，在相同情形下计算512个特征点的描述子时，SURF耗时335ms,BRIEF仅8.18ms；匹配SURF描述子需28.3ms，BRIEF仅需2.19ms。在要求不太高的情形下，BRIEF描述子更容易做到实时。

改进BRIEF算法—rBRIEF（Rotation-AwareBrief）

（1）steered BRIEF（旋转不变性改进）

在使用oFast算法计算出的特征点中包括了特征点的方向角度。假设原始的BRIEF算法在特征点SxS（一般S取31）邻域内选取n对点集。经过旋转角度θ旋转，得到新的点对，在新的点集位置上比较点对的大小形成二进制串的描述符。这里需要注意的是，在使用oFast算法是在不同的尺度上提取的特征点。因此，在使用BRIEF特征描述时，要将图像转换到相应的尺度图像上，然后在尺度图像上的特征点处取SxS邻域，然后选择点对并旋转，得到二进制串描述符。

（2）rBRIEF-改进特征点描述子的相关性

使用steeredBRIEF方法得到的特征描述子具有旋转不变性，但是却在另外一个性质上不如原始的BRIEF算法。是什么性质呢，是描述符的可区分性，或者说是相关性。这个性质对特征匹配的好坏影响非常大。描述子是特征点性质的描述。描述子表达了特征点不同于其他特征点的区别。我们计算的描述子要尽量的表达特征点的独特性。如果不同特征点的描述子的可区分性比较差，匹配时不容易找到对应的匹配点，引起误匹配。

为了解决描述子的可区分性和相关性的问题，ORB使用统计学习的方法来重新选择点对集合。

首先建立300k个特征点测试集。对于测试集中的每个点，考虑其31x31邻域。这里不同于原始BRIEF算法的地方是，这里在对图像进行高斯平滑之后，使用邻域中的某个点的5x5邻域灰度平均值来代替某个点对的值，进而比较点对的大小。这样特征值更加具备抗噪性。另外可以使用积分图像加快求取5x5邻域灰度平均值的速度。

从上面可知，在31x31的邻域内共有(31-5+1)x(31-5+1)=729个这样的子窗口，那么取点对的方法共有M=265356种，我们就要在这Ｍ种方法中选取256种取法，选择的原则是这256种取法之间的相关性最小。怎么选取呢？

1. 在300k特征点的每个31x31邻域内按M种方法取点对，比较点对大小，形成一个300kxM的二进制矩阵Q。矩阵的每一列代表300k个点按某种取法得到的二进制数。
2. 对Q矩阵的每一列求取平均值，按照平均值到0.5的距离大小重新对Q矩阵的列向量排序，形成矩阵T。
3. 将T的第一列向量放到R中。
4. 取T的下一列向量和R中的所有列向量计算相关性，如果相关系数小于设定的阈值，则将T中的该列向量移至R中。
5. 按照上一步的方式不断进行操作，直到R中的向量数量为256。

这就是rBRIEF算法。