SIFT角点学习

SIFT(Scale-invariant feature transform)

本文为笔记原理详见资料1
资料1
资料2
资料3

高斯模糊

	SIFT算法是在不同的尺度空间上查找关键点，而尺度空间的获取需要使用高斯模糊来实现，Lindeberg等人已证明高斯卷积核是实现尺度变换的唯一变换核，并且是唯一的线性核。

1. 二维高斯函数
   $\sigma :正太分布标准差，值越大图像越模糊$
   $r:模糊半径：元素到模板中心的距离$
   $$G(r)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-r^2/2{\sigma }^2}$$
   如果二维模板大小为m×n，则模板上的元素(x,y)对应的高斯计算公式为：
   $$G(x,y)=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(x-m/2{)}^2+(y-n/2{)}^2}{2{\sigma }^2}}$$
   其效果如下：
   

高斯金字塔的构建

尺度空间在实现时使用高斯金字塔表示，高斯金字塔的构建分为两部分：
（1）：对图像做不同尺度的高斯模糊
（2）：对图像做降采样（缩放）,将原始图像不断降阶采样，得到一系列大小不一的图像，由大到小，从下到上构成的塔状模型。原图像为金子塔的第一层，每次降采样所得到的新图像为金字塔的一层(每层一张图像)，每个金字塔共n层。金字塔的层数根据图像的原始大小和塔顶图像的大小共同决定，其计算公式如下：
n = l o g 2 { m i n ( M , N ) } − l o g 2 t ( 金 字 塔 共 n 层 , t 为 最 小 维 数 图 像 的 对 数 值 ) n=log_2\{min(M,N)\}-log_2t(金字塔共n层,t为最小维数图像的对数值) n=log2​{min(M,N)}−log2​t(金字塔共n层,t为最小维数图像的对数值)
高斯金字塔在简单降采样的基础上加上了高斯滤波，将图像金字塔每层的一张图像使用不同参数做高斯模糊，使得金字塔的每层含有多张高斯模糊图像，将金字塔每层多张图像合称为一组(Octave)，金字塔每层只有一组图像，组数和金字塔层数n相等，，每组含有多张(也叫层Interval)图像。

尺度

2. 尺度空间理论
   模拟人在距离目标由近到远时目标在视网膜上的形成过程，在图像信息处理模型中引入一个被视为尺度的参数，通过连续变化尺度参数获得多尺度下的尺度空间表示序列，对这些序列进行尺度空间主轮廓的提取，并以该主轮廓作为一种特征向量，实现边缘、角点检测和不同分辨率上的特征提取等。
   （1）：尺度空间的表示
   $L(x,y,\sigma ):表示卷积$
   $G(x,y,\sigma )=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(x-m/2{)}^2+(y-n/2{)}^2}{2{\sigma }^2}}$
   $$L(x,y,\sigma )=G(x,y,\sigma )\ast I(x,y)$$

   (1)尺度空间
   $O:金字塔的的塔数（Y方向）$
   $w,h:图像的宽和高$
   $$O=[lo{g}_2(min(w,h))]-3$$
   $S每层塔的组数（X方向），一般选3～5$：
   这样由$(O,S)$就构成了尺度空间.
   (2)如何确定每组中每个图像对应的的尺度？($\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(x-m/2{)}^2+(y-n/2{)}^2}{2{\sigma }^2}}$中的$\sigma$):
   令$k=2^{\frac{1}{S}}$:
   对于同一组内基准层尺度与各尺度关系：
   $i表示第i层$
   $$\sigma (i)={\sigma }_0{2}^{\frac{i}{S}}$$
   对于不同组（不同的塔层）内基准层尺度与各尺度关系:
   $$\sigma (O,S,i)={\sigma }_0{2}^{O+\frac{i}{S}}$$

关键点检测

1. 高斯差分金字塔（DOG空间）

   2002年Mikolajczyk在详细的实验比较中发现尺度归一化的高斯拉普拉斯函数的极大值和极小值同其它的特征提取函数。在实际计算时，使用高斯金字塔每组中相邻上下两层图像相减，得到高斯差分图像，如图3.3所示，进行极值检测。

其中拉普拉斯函数和高斯差分函数具有如下关系：
${\sigma }^2{▽}^{2}G:利亚谱诺夫函数$
$D(x,y,\sigma ):高斯差分函数$
$$\begin{gathered} {\sigma }^2{▽}^{2}G=\frac{\partial G}{\partial \sigma }\approx \frac{G(x,y,k\sigma )-G(x,y,\sigma )}{k\sigma -\sigma } \\ G(x,y,k\sigma )-G(x,y,\sigma )\approx (k-1){\sigma }^2{▽}^{2}G \end{gathered}$$
求极值时使用高斯差分算子带起拉普利斯算子进行极值检测，如下:
$\sigma (i)={\sigma }_0\sqrt{2^{i/S}-2^{(i-1)/S}},i\in [0,..,s+2]:尺度空间坐标$
$$\begin{gathered} D(x,y,\sigma )=(G(x,y,k\sigma )-G(x,y,\sigma ))\ast I(x,y) \\ =L(x,y,k\sigma (i+1))-L(x,y,\sigma (i)) \end{gathered}$$

2. 空间极值点检测
   
   为了寻找DoG函数的极值点，每一个像素点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。如图3.4所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。
   由于要在相邻尺度进行比较，如图3.3右侧每组含4层的高斯差分金子塔，只能在中间两层中进行两个尺度的极值点检测，其它尺度则只能在不同组中进行。为了在每组中检测S个尺度的极值点，则DOG金字塔每组需S+2层图像，而DOG金字塔由高斯金字塔相邻两层相减得到，则高斯金字塔每组需S+3层图像，实际计算时S在3到5之间。
   尺度详解

3. 关键点定位
   通过拟合三维二次函数来精确确定关键点的位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点(因为DoG算子会产生较强的边缘响应)，以增强匹配稳定性、提高抗噪声能力。
   （1）： 离散空间的极值点并不是真正的极值点，图4.1显示了二维函数离散空间得到的极值点与连续空间极值点的差别。利用已知的离散空间点插值得到的连续空间极值点的方法叫做子像素插值。
   
   为了提高关键点的稳定性，需要对尺度空间DoG函数进行曲线拟合。利用DoG函数在尺度空间的Taylor展开式(拟合函数)为：
   $X=(x,y,\sigma {)}^T$
   $$D(X)=D+\frac{\partial D^T}{\partial X}X+\frac{1}{2}{X}^T\frac{{\partial }^{2}D}{\partial X^2}X=0$$
   解得极值点：
   $$\hat{X}=-\frac{{\partial }^2{D}^{-1}}{\partial X^2}\frac{\partial D}{\partial X}$$
   极值为：
   $$D(\hat{X})=D+\frac{1}{2}\frac{\partial D^T}{\partial X}\hat{X}$$
   $X=(x,y,\sigma {)}^T$代表相对插值中心的偏移量，当它在任一维度上的偏移量大于0.5时（即x或y或$\sigma$），意味着插值中心已经偏移到它的邻近点上，所以必须改变当前关键点的位置。同时在新的位置上反复插值直到收敛；也有可能超出所设定的迭代次数或者超出图像边界的范围，此时这样的点应该删除，一般迭代5次。
   $|D(x)|$的值:过小容易受噪声影响，变得不稳定，需要对$|D(x)|$进行限幅一般为0.03或者0.04.
   （2）：消除边缘效应
   DOG算子会产生较强的边缘响应，需要剔除不稳定的边缘响应点。获取特征点处的Hessian矩阵，主曲率通过一个2x2 的Hessian矩阵H求出：
   $$H=\left[\begin{array}{cc}{D}_{xx}& D_{xy}\\ D_{xy}& D_{yy}\end{array}\right]$$
   $Det(H)=D_{xx}{D}_{yy}-(D_{xy})=\alpha +\beta$
   $tr(H)=D_{xx}-D_{yy}=\alpha \beta$
   令$\alpha =r\beta$求：$\frac{tr(H)}{Det(H)}$
   $$\frac{tr(H)}{Det(H)}=\frac{(r+1{)}^2}{r}$$
   D的主曲率和H的特征值成正比，令为α最大特征值，β为最小的特征值，则公式的值在两个特征值相等时最小，随着的增大而增大。值越大，说明两个特征值的比值越大，即在某一个方向的梯度值越大，而在另一个方向的梯度值越小，而边缘恰恰就是这种情况。所以为了剔除边缘响应点，需要让该比值小于一定的阈值，因此，为了检测主曲率是否在某域值r下，只需检测$\frac{tr(H)}{Det(H)}<\frac{(r+1{)}^2}{r}$成立时将关键点保留，反之删除。

关键点方向检测

为了使描述符具有旋转不变性，需要利用图像的局部特征为给每一个关键点分配一个基准方向。
在DOG金字塔中检测出关键点之后,采集关键点邻域内（4×4区域8个方向的特征向量）的灰度变化，使算子具有尺度和方向不变性。

1. 确定计算描述子所需的图像区域。
   对梯度的求取应在特征点对应的高斯图像上进行。
   梯度的模值m(x,y)按的高斯分布加成，按尺度采样的3σ原则。
   邻域窗口半径为：$1\times 1.5\times {\sigma }_{oct}(i)({\sigma }_{oct}={\sigma }_0{2}^{i/S}(i\in (1,..,S+2)))$。
   梯度模值为：$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+1,y)-L(x-1,y)}{L(x,y+1)-L(x,y-1)})$
2. 使用直方图统计邻域内像素的梯度和方向
   梯度直方图将0~360度的方向范围分为36个柱(bins)，其中每柱10度。如图5.1所示，直方图的峰值方向代表了关键点的主方向，(为简化，图中只画了八个方向的直方图)。
   方向直方图的峰值则代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向作为该关键点的辅方向。因此，对于同一梯度值的多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被创建但方向不同。仅有15％的关键点被赋予多个方向，但可以明显的提高关键点匹配的稳定性。实际编程实现中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点，并且，离散的梯度方向直方图要进行插值拟合处理，来求得更精确的方向角度值，检测结果如图5.2所示。
   

关键点特征描述

通过以上步骤，对于每一个关键点，拥有三个信息：位置、尺度以及方向。接下来就是为每个关键点建立一个描述符，用一组向量将这个关键点描述出来，使其不随各种变化而改变，比如光照变化、视角变化等等。这个描述子不但包括关键点，也包含关键点周围对其有贡献的像素点，并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。
描述子使用在关键点尺度空间内$4\times 4$的窗口中计算的8个方向的梯度信息，共$4\ast 4\ast 8=128$维向量表征。表示步骤如下：

1. 确定计算描述子所需的图像区域
   对梯度的求取应在特征点对应的高斯图像上进行。将关键点附近的邻域划分为d*d(Lowe建议d=4)个子区域，每个子区域做为一个种子点，每个种子点有8个方向。每个子区域的大小与关键点方向分配时相同，即每个区域有个$3{\sigma }_{oct}$子像素，为每个子区域分配边长为$3{\sigma }_{oct}$的矩形区域进行采样(为了简化计算取其边长为$3{\sigma }_{oct}$，并且采样点宜多不宜少)。考虑到实际计算时，需要采用双线性插值，所需图像窗口边长为。在考虑到旋转因素(方便下一步将坐标轴旋转到关键点的方向)，如下图6.1所示，实际计算所需的图像区域半径为：
   $$r=\frac{3{\sigma }_{oct}\times \sqrt{2}\times (d+1)}{2}(四舍五入)$$
   

2. 将坐标轴旋转为关键点的方向，以确保旋转不变性，如6.2所示。
   
   $$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]=\left[\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right]=\left[\begin{array}{c}x\\ y\end{array}\right](x,y\in [-r,r])$$

3. 将邻域内的采样点分配到对应的子区域内.
   将子区域内的梯度值分配到8个方向上，计算其权值。旋转后的采样点坐标在半径为r的圆内被分配到的$d\times d$子区域，计算影响子区域的采样点的梯度和方向，分配到8个方向上。旋转后的采样点$(x^{\prime },y^{\prime })$落在子区域的下标为:
   $$\left[\begin{array}{c}{x}^{\prime \prime }\\ y^{\prime \prime }\end{array}\right]=\frac{1}{3{\sigma }_{oct}}\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right]+\frac{d}{2}$$
   子区域的像素的梯度大小按$\sigma =0.5d$的高斯加权计算，即:
   $$w=m(a+x,b+y)e^{-\frac{(x^{\prime }{)}^2+(y^{\prime }{)}^2}{2\times (0.5d{)}^2}}$$
   其中a，b为关键点在高斯金字塔图像中的位置坐标。

4. 采用双线性插值的方法插出网格点的值：
   参考
   
   将由式(6-3)所得采样点在子区域中的下标$(x^{\prime \prime },y^{\prime \prime })$(图中蓝色窗口内红色点)线性插值，计算其对每个种子点的贡献。如图中的红色点，落在第0行和第1行之间，对这两行都有贡献。对第0行第3列种子点的贡献因子为$dr$，对第1行第3列的贡献因子为$1-dr$，同理，对邻近两列的贡献因子为dc和1-dc，对邻近两个方向的贡献因子为$do$和$1-do$。则最终累加在每个方向上的梯度大小为：
   $k,m,n=0,1根据实际位置选取$
   $w:红色点的梯度值$
   $$weight=w\ast d{r}^k\ast (1-dr{)}^{1-k}\ast d{c}^m\ast (1-dr{)}^{1-m}d{o}^n\ast (1-do{)}^{1-o}$$

5. 为了去除光照变化的影响
   如上统计的448=128个梯度信息即为该关键点的特征向量。为了去除光照变化的影响，需要对它们进行归一化处理，对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到，所以也能去除。得到的描述子向量为，归一化后的特征向量为则如上统计的448=128个梯度信息即为该关键点的特征向量。特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到，所以也能去除。得到的描述子向量为$H=(h_1,h_2,...,h_{128})$，归一化后的特征向量$L=(l_1,l_2,...,l_{1}28)$为则
   $$l_i=\frac{h_i}{\sqrt{{\sum }_{j=1}^{128}}},j=1,2,3,..$$

6. 描述子向量门限
   非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此设置门限值(向量归一化后，一般取0.2)截断较大的梯度值。然后，再进行一次归一化处理，提高特征的鉴别性。

7. 按特征点的尺度对特征描述向量进行排序
   资料
   为了排除因为图像遮挡和背景混乱而产生的无匹配关系的关键点，SIFT的作者Lowe提出了比较最近邻距离与次近邻距离的SIFT匹配方式：取一幅图像中的一个SIFT关键点，并找出其与另一幅图像中欧式距离最近的前两个关键点，在这两个关键点中，如果最近的距离除以次近的距离得到的比率ratio少于某个阈值T，则接受这一对匹配点。因为对于错误匹配，由于特征空间的高维性，相似的距离可能有大量其他的错误匹配，从而它的ratio值比较高。显然降低这个比例阈值T，SIFT匹配点数目会减少，但更加稳定，反之亦然。
   　　Lowe推荐ratio的阈值为0.8，但作者对大量任意存在尺度、旋转和亮度变化的两幅图片进行匹配，结果表明ratio取值在0. 4~0. 6 之间最佳，小于0. 4的很少有匹配点，大于0. 6的则存在大量错误匹配点，所以建议ratio的取值原则如下:
   　　ratio=0. 4：对于准确度要求高的匹配；
   　　ratio=0. 6：对于匹配点数目要求比较多的匹配；
   　　ratio=0. 5：一般情况下