图像配准系列之Sift特征点检测与匹配(1)
在讲解Sift特征点的检测与匹配之前，先讲一下本人对图像配准与特征点的理解。
1. 图像配准
图像配准，就是找到一幅图像与另一幅图像中相同位置点的一一映射关系，然后根据点映射关系，对其中的一幅图像（通常称为浮动图像或待配准图像）进行空间坐标变换（也称为像素重采样），使其与另一幅图像（通常称为参考图像）位置匹配。
如下图所示，浮动图像与参考图像中具有相似的蓝色区域与黑点，但两者的位置并不匹配，也即相同部位在两张图像上的坐标位置并不一样，浮动图像相对于参考图像来说具有旋转与平移的位置偏差，图像配准的目的就是为了矫正浮动图像的位置，消除这些位置偏差，从而使浮动图像与参考图像位置匹配。
基于Sift特征点匹配的图像配准流程如下图：
2. 特征点
顾名思义，特征点是“点+特征”的组合，“点”代表了特征点的位置，通常为图像中灰度值相对其周围区域发生剧烈变化的点，比如纹理剧烈变化的点、拐角点、直线与直线的交点，以及单纯区域中的孤立点等。“特征”则对该点位置的特性进行描述，其不仅包含了该点本身的特性，还包含了该点周围邻域点的特性，“特征”通常使用一个n维向量来表征。
比如Sift特征点，每个Sift特征点包含了灰度值相对其周围区域发生剧烈变化的点位置，以及用于描述该点与其邻域特征的一个128维向量。
3. Sift特征点检测匹配的原理与过程
Sift特征点的检测步骤大致分为以下几个步骤：
(1) 建立高斯金字塔图像
Sift算法使用的高斯金字塔图像的建立过程，通常分为降采样和高斯卷积两个步骤。高斯金字塔分为多个大层，每个大层又分为多个小层，不同大层包含图像的尺寸与高斯卷积参数都不一样，然而属于同一个大层的不同小层图像，其尺寸是一样的，但高斯卷积参数不一样。
自底向上，第0大层的宽和高与原图一致，第1大层的宽、高分别是第0大层宽、高的1/2，第2大层的宽、高分别是第1大层宽、高的1/2，依此类推......
a. 降采样，也就是对图像取奇数行、列得到的向下采样图像。
> 第0大层的图像与原始图像尺寸一样，不需要降采样。
> 第1大层的第0小层图像由第0大层的倒数第3小层图像（上图中的倒数第3小层为自底向上数的第1小层）降采样得到，不需要降采样。第1大层的其它小层的图像则由它的第0小层图像经过高斯卷积得到。
> 第2大层的第0小层图像由第1大层的倒数第3小层图像降采样得到，不需要降采样。第2大层的其它小层的图像则由它的第0小层图像经过高斯卷积得到。
> 第3大层的第0小层图像由第2大层的倒数第3小层图像降采样得到，不需要降采样。第3大层的其它小层的图像则由它的第0小层图像经过高斯卷积得到。
依此类推......
b. 高斯卷积，使用高斯卷积核对图像进行高斯模糊。通常使用5*5的高斯卷积核，5*5的卷积核的计算式如下，其中σ为标准差，σ越大，权重分布越均匀，滤波效果越好，同时图像越模糊，反之σ越小，权重分布越偏向于窗口中心点，滤波效果越差，同时图像越能保留其原有清晰度。
假设高斯金字塔图像有M个大层，每个大层有N个小层，对于每一小层，其所属的大层序列号为m(0≤m≤M-1)，小层序列号为n(0≤n≤N-1)，初始参数为σ0(通常取1.6)，那么每一小层的高斯卷积参数σ可按下式计算(注意：每个大层的第0小层都不需要做高斯卷积)：
上式中为什么x=N-3，在下文我们再详细说明。
(2) 计算高斯金字塔差分图像
高斯金字塔差分图像，又称为DOG金字塔图像，是由以上第(1)步所得的高斯金字塔图像计算得来的，在每一个大层中，计算所有相邻小层的图像差值，即为差分图像，所得的所有差分图像组成高斯金字塔差分图像。如下图所示：
(3) 检测极值点
极值点，直观来说就是上文说到的灰度值相对其周围区域发生剧烈变化的点。使用上述第(2)步得到的高斯金字塔差分图像来判断一个点是否是极值点（高斯金字塔差分图像同样具有大层与小层的概念），比如在下图中，要判断点A（红点）是否是极值点，则不仅取点A的周围邻域8个点与点A的像素值进行比较，还取点A所在小层的上、下两个相邻小层中对应位置的3*3点进行比较，因此总共需要取8+3*3+3*3=26个点与点A的像素值进行比较，如果点A的像素值在26个点中是最大值或者最小值，则判定点A是极值点。
上文已经提到，在高斯金字塔图像中，假设每个大层有N个小层，那么由N个小层图像计算所得的差分图像个数为N-1。然而判断一个点是否是极值点，不仅取差分图像上该点的邻域点进行比较，还取该点所在小层的上、下两个相邻小层的差分图像中对应位置的3*3点进行比较，因此在每个大层中最底部小层与最顶部小层是无法检测极值点的。由此可知，对于高斯金字塔差分图像的每一个大层，需要检测极值点的小层个数为N-1-1-1=N-3，x=N-3就是这么来的。
经过降采样和高斯模糊之后，得到的高斯金字塔图像是离散的，同样得到的高斯金字塔差分图像也是离散的，导致以上求得的极值点很可能存在位置偏差，因此需要对极值点位置作进一步精确定位。假设检测到的极值点为点X(x,y,σ)，修正的极值点为点X0(x0,y0,σ0)，点X的像素值为f(X)，将f(X)在点X0处作三元二阶泰勒展开并忽略高阶项，有：
上式(1)中的h(△X)为：
对h(△X)中的△X求导，得到：
又因为H是对称矩阵，H的转置与H是完全相等的，因此有：
对以上(1)式中的△X求导，并把(3)式代入得到：
这里需要注意，向量对向量的求导遵从以下规则：
A和x都是向量，A与x的矩阵相乘对x的导数为A的转置：
因此有：
为了求极值，令以上(4)式等于0，那么解得：
将以上(5)式代入(1)式，得到：
在以上的(5)式中，△X为已检测到极值点X与修正极值点X0的坐标偏移量，根据(5)式，我们只要求得点X(x,y,σ)的所有一阶偏导数与二阶偏导数，就可以求得坐标偏移量△X。通常使用差分梯度来近似求得偏导数的值。假设检测到的极值点X所在的小层为图像Icur，其上一小层为Inext，下一小层为Ipre，那么点X的偏导数计算如下式（同一小层中所有像素点的σ是一样的，因此对σ的偏导数表现为使用不同小层的点来计算差分）：
根据上式求解得△X(△x1,△x2,△x3)的值之后，判断如果是否△x1或△x2或△x3的值超出一定范围，如果超出一定范围则舍弃该解。如果没有超出一定范围，则使用△X来修正X的位置，也即求得点X0的位置：
修正X的位置是一个迭代的过程：当求得X0的位置之后，X0的位置比X更精准，此时把点X0当作点X，重复(5)式的计算得到新的△X并修正X位置。通常迭代次数为5次，超出5次则退出迭代。
此外，如果根据以上(6)式求得的f(X)值小于一定阈值（阈值通常取经验值0.03或0.04/N，N为每个大层的小层数），则认为该点不稳定，也将其剔除。
首先将f(X)归一化到区间[0,1]，然后再比较是否：
高斯差分算子会产生较强的边缘响应，极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。为了进一步确保极值点的稳定，接下来还要剔除具有较强边缘响应的点。经过以上步骤，已经得到比较精确的极值点X，我们通过以下步骤来剔除边缘响应点：
首先，计算点X的海塞矩阵：
假设α和β分别为H的特征值，且α为较大特征值，并设α=γβ，接着计算H的对角元素和Tr(H)与行列式值Det(H)：
如果Det(H)<0则舍弃X。
然后判断下式是否满足，如果不满足也舍弃X。下式的γ0通常取经验值10。
(4) 计算描述子
① 获取极值点对应的高斯图像
经过上述步骤，我们得到了精确稳定的极值点，也即得到了若干极值点的坐标X(x,y,σ)，坐标中的σ表示其尺度，σ表明了该极值点属于DOG金字塔图像中的哪一大层哪一小层，从而确定了该极值点对应高斯金字塔图像中的哪一大层哪一小层图像Lσ。高斯金字塔图像与DOG金字塔图像的层对应关系如下图：
这时我们需要在Lσ上统计以X(x,y,σ)为中心，以3*1.5σ为半径的圆内所有点的梯度大小与梯度方向，如下图所示：
圆形区域内每个点的梯度大小g(x,y)与梯度方向θ(x,y)的计算如下式：
为了防御噪声的干扰，计算得到圆形区域内所有点的梯度大小与方向之后，以极值点为中心，以1.5σ为高斯分布的标准差，对圆形区域内所有点的梯度大小进行高斯加权（越靠近中心点权重越大）：
上式中，高斯加权的权重需除以区域内所有点权重之和进行归一化，d为圆形区域内每个邻域点与极值点的像素距离。
② 统计梯度直方图
θ的取值范围是0°~360°，我们将0°~360°平均分成36个区间，统计圆形区域内每个点的梯度方向属于哪个区间，然后把该点的梯度大小累加到该区间的幅度（也即每个区间的幅度等于该区间内所有点的梯度大小之和），直观来看就是一个统计直方图：
为进一步防御噪声的干扰，还需要对梯度直方图的幅值进行平滑处理。假设从区间0°~9°到区间350°~360°，以上步骤所得的梯度直方图赋值分别为h1~h36，那么平滑公式为：
根据以上公式平滑直方图，有可能下表超出1~36的范围，这时可以把1~36的幅值看作一个圆形来循环取值：
由此，我们得到经过平滑之后的36个直方图幅值H1~H36，分别对应0°~9°、10°~19°、20°~29°、......、350°~360°这36个角度区间，直方图幅值最大（直方图峰值）的角度区间定为该极值点的主方向，其它所有直方图幅值大于“直方图峰值*80%”的角度区间，定为辅方向，因此一个极值点有可能产生多个特征点：“一个主方向的特征点+几个辅方向的特征点”。对应同一个极值点的多个特征点，除了方向不一样，位置和尺度都是一样的。
③ 梯度直方图的抛物线插值
很容易可以发现，上一步得到的主方向与辅方向，都是某一个角度区间，而不是某一个具体确定的角度值，所以我们还得由角度区间进一步确定具体的角度值。Sift使用抛物线插值的方法来确定角度值。仅对满足以下两个条件的角度区间进行抛物线插值，对于不完全满足这两个条件的角度区间（方向）则将其舍弃：
a. 是主方向或辅方向；
b. 该区间的梯度幅值大于其两侧区间的梯度幅值。
假设满足以上条件的区间序列号为k(2≤k≤35)，其左侧区间序列号为k-1，右侧区间序列号为k+1（k-1、k、k+1均对应角度区间的中点），如下图所示：
于是我们得到平面坐标系上的3个点(k-1，Hk-1)，(k，Hk)，(k+1，Hk+1)。假设抛物线的解析式为：
将以上3个坐标点代入解析式解得：
由抛物线的特性可知，当x=-b/(2a)时，抛物线取得最大值或最小值，因此抛物线顶点对应的k值为：
上式的Δk为k的偏移量，k+Δk就是我们要求的抛物线插值，我们知道：
因此可以轻易证明：
也就是说，kmax=k+Δk仍然属于第k个角度区间。假设第k个角度区间的角度范围是α0~α1，k值对应的是该区间的中点α0+(α1-α0)/2=(α0+α1)/2，那么kmax对应的角度α为：
上式中α就是我们要确定的主方向或辅方向的角度值。
接下来还有描述子（128维向量）的生成、匹配，误匹配的剔除、计算变换矩阵，像素重采样等步骤，这个篇幅有点长了，本文就写到这里，我们将在下一篇文章中继续讲解接下来的内容。
微信公众号如下，欢迎扫码关注，欢迎私信技术交流：
