SIFT是由Lowe[2004]开发的一个算法，用于提取图像中的不变特征。称它为变换的原因是，它会将图像数据变换为相对于局部图像特征的尺度不变坐标。SIFT是本章到目前为止讨论的最复杂的特征检测和描述方法。本节中用到了大量由实验确定的参数。因此，与前面介绍的许多方法不同，SIFT具有很强的探索性，因为我们目前所具备的知识还不能让我们将各种方法组合为一个能够求解单重方法所不能求解的问题的“系统”。因此，我们不得不通过实验确定控制更复杂系统性能的各种参数的相互作用。SIFT特征（称为关键点）对图像尺度和旋转是不变的，并且对仿射失真、三维视点变化、噪声和光照变化具有很强的鲁棒性。SIFT的输入是一幅图像，输出是一个n维特征向量，向量的元素是不变的特征描述子。
Lowe将SIFT算法分为如下四步：
（1）尺度空间极值检测：搜索所有尺度上的图像位置。通过高斯差分函数来识别潜在的对于尺度和旋转不变的关键点。
（2）关键点定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。
（3）关键点方向确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进项变换，从而保证了对于这些变换的不变性。
（4）关键点描述：在每个关键点周围的邻域内，在选定的尺度上测量图像局部的梯度。这些梯度作为关键点的描述符，它允许比较大的局部形状的变形或光照变化。
1.1.1尺度空间
SIFT算法的第一阶段是找到对尺度变化不变的图像位置。将图像利用不同的尺度参数进行平滑后得到一堆图像，目的是模拟图像的尺度减小出现的细节损失。在SIFT中，用于实现平滑的是高斯核，因此尺度参数是标准差。标准差的大小决定图像的平滑程度，大尺度对应图像的概貌特征，小尺度对应图像的细节特征。大的值对应粗糙尺度(低分辨率)，小的值对应精细尺度(高分辨率)。因此灰度图像f(x,y)的尺度空间L(x,y,σ)是f与一个可变尺度高斯核G(x,y,σ)的卷积：
式中，尺度由参数σ控制，G的形式如下：
输入图像f(x,y)依次与标准差为σ，kσ，k2σ，k3σ，...，的高斯核卷积，生成一堆由一个常量因子k分隔的高斯滤波（平滑）图像。
Q1:怎么理解平滑建立的尺度空间？
对于一副图像，近距离观察和远距离观察看到的图像效果是不同的，前者比较清晰，通过前者能看到图像的一些细节信息，后者比较模糊，通过后者能看到图像的一些轮廓的信息，这就是图像的尺度，图像的尺度是自然存在的，并不是人为创造的。尺度空间中各尺度图像的模糊程度逐渐变大，能够模拟人在距离目标由近到远时目标在视网膜上的形成过程。在平滑后图像上的一点，包含的是原图中一块像素的信息，属于语义层面的多尺度。
1.1.2高斯金字塔
图像金字塔广泛应用于各种视觉应用中。图像金字塔是图像的集合，它是由原始图像产生，连续降采样，直到达到一些期望的停止点。常出现的两种金字塔是：高斯金字塔和拉普拉斯金字塔。高斯金字塔用于降采样，当我们要从金字塔中较低的图像重构上采样图像时，需要拉普拉斯金字塔。SIFT中的高斯金字塔与最原始的高斯金字塔有些区别，因为它在构造尺度空间时，将这些不同尺度图像分为了O个Octave、每个Octave又分为了S层。
上面一行原始信号类似于一条线的灰度剖面，下面是在相同尺度σ下的LoG结果。可以看出总的LoG曲线其实是两条边界上LoG函数的结果的叠加，当两条边界足够接近时，可以将这个原始信号看作一个blob，两个LoG函数的结果也合二为一得到一个极值点，这个极值点对应blob的中心。所以边缘检测对应的是LoG的过零点，而斑点检测对应的是LoG的极值点。
检测blob时，需要选择合适的尺度，然后LOG可以通过检测极值点来检测blob。上图第一行为原信号对不同尺度LoG的响应，可以发现随着尺度的不断增大，LoG曲线由双波谷逐渐融合成单波谷，但是响应的幅值越来越弱。这是因为，随着尺度的增大，LoG算子的最大幅度逐渐减小，导致响应也随着尺度的增大而减小。因此应该使用σ2对LoG进行归一化，即σ2▽2G。上图第二行是归一化后的拉普拉斯响应。归一化后选择产生最强响应的尺度，在该尺度上对应的极值就是blob的中心位置。理论表明，对于一个圆形blob，当二维LoG算子的零点值曲线和目标圆形边缘重合时取得最强响应。
即令LoG=0，
得到σ=r/√2。使用LoG寻找斑点时不仅在图像上寻找极值点，还要求在尺度空间上也是极值点。也就是检测在尺度空间和图像空间都是极值的点，就是blob区域的中心点。
1.1.4高斯差分(DOG)金字塔
而LoG函数为
两者仅相差一个σ，即
根据极限又可以得到
因此，当k≈1时，有
两个不同方差的高斯函数相减记作DoG，和尺度归一化的LoG只相差一个倍数k-1，这个倍数是一个常数，不影响极值点检测。DoG计算速度快，LoG更精确。因此，我们用高斯差DoG近似高斯拉普拉斯LoG，具体做法就是将相邻两层的高斯尺度图像相减：
OpenCV构建DoG金字塔源码
并行计算，共nOctaves个Octave，每个Octave有nOctaveLayers+3层，共nOctaves*(nOctaveLayers+2)个差分计算。
由高斯金字塔得到每组n+3层图像，所有相邻两幅高斯滤波图像得到n+2个差函数D(x,y,σ)。将这些差函数视为图像，图像的细节水平随着我们上调尺度空间而下降。下图显示了SIFT查找图像D(x,y,σ)中的极值的过程。
在D(x,y,σ)图像中的每个位置（显示为黑色x），将该位置的像素值与其在当前图像中的8个相邻像素值及其在上方和下方图像中的9个相邻像素值，即26个相邻像素，进行比较。如果该位置的值大于其所有相邻像素的值，或小于所有相邻像素的值，那么该位置被选为极值（最大值或最小值点）。在一个Octave的第一个（最后一个）尺度中检测不到任何极值，因为它没有相同大小的下方（上方）尺度图像。
OpenCV查找初始关键点源码
在DoG尺度空间中找初始极值点，是通过比较每一个像素与其26个邻域像素相比较。如果该位置的值大于其所有相邻像素的值，或小于所有相邻像素的值，那么该位置被选为极值。OpenCV中对应的是findScaleSpaceExtremaT类，其中包括了极值计算，调整极值，梯度计算等。下面为极值计算部分。
一个连续函数被取样时，它真正的最大值或最小值实际上可能位于样本点之间。得到接近真实极值（亚像素精度）的一种方法时，首先再数字函数中的每个极值点处拟合一个内插函数，然后在内插后的函数中查找改进精度后的极值位置。SIFT用D(x,y,σ)的泰勒级数展开的线性项和二次项，把原点移至被检测的样本点。这个公式的向量形式为
式中，D及其导数是在这个样本点处计算的，x=(x,y,σ)T是这个样本的偏移量，▽是我们熟悉的梯度算子，
H是海森矩阵，
取D(x)关于x的导数并令其为零，可求得极值的位置x，即
SIFT使用极值位置函数值D(x)来剔除具有低对比度的不稳定极值，其中D(x)为
Lowe的实验结果称，若所有图像的值都在区间[0,1]内，则D(x)小于0.03的任何极值都会被剔除。这剔除了具有低对比度和/或局部化较差的关键点。
1.2.3消除边缘响应
H的特征值与D的曲折度成正比。如对哈里斯-斯蒂芬斯(HS)角检测器的解释那样，我们可以避免直接计算特征值，H的迹等于特征值之和，H的行列式等于特征值之积。设α和β分别是H的最大特征值和最小特征值，有
注意，H是对称的，且大小为2×2。若行列式是负的，则不同曲折度具有不同的符号，且讨论的关键点不可能是一个极值，因此要丢弃它。令r表示最大特征值和最小特征值之比，即α=rβ，则
当特征值相等时，上式出现最小值，当r大于1时，随着r的增大而增大。因此，要检查低于某个阈值的r的主曲折度之比，只需检查
这个计算很简单。Lowe报告的实验结果中使用了值r=10，这意味着消除了曲折度之比大于10的关键点。
OpenCV源码adjustLocalExtrema
包括计算亚像素精度极值+剔除低对比度的不稳定极值+剔除曲折度较大的点，基于Lowe论文的Section4
OpenCV源码calcOrientationHist
1.校正旋转主方向，确保旋转不变性。为了保证特征矢量的旋转不变性，要以特征点为中心，在附近邻域内将坐标轴旋转θθ（特征点的主方向）角度，即将坐标轴旋转为特征点的主方向。旋转后邻域内像素的新坐标为：
2.生成描述子，最终形成一个128维的特征向量。
以关键点为中心取16×16像素的一个区域，并使用像素差计算这个区域中每个点处的梯度幅度和方向。然后使用标准差等于这个区域大小一般的高斯加权函数来分配一个权重，将这个权重乘以每个点处的梯度幅度。高斯加权函数在图中显示为一个圆，权重随着到中心的距离增大而减小。这个函数的目的是，减少位置变化很小时描述子的突然变化。
将16×16的区域分为16个4×4区域，每个4×4区域中有16个方向，将所有的梯度方向量化为8个相隔45°的方向。SIFT不是把一个方向值分配给其最接近的容器，而是进行内插运算，根据这个值到每个容器中心的距离，按比例地在所有容器中分配直方图的输入。例如，容器的中心为[22.5,67.5,112.5,...,337.5]，每个输入需要乘以权重1-d，d是从这个值到容器中心的最短距离，度量单位是直方图间隔，最大的可能距离是1。假设某个方向值是22.5°，到第一个容器中心的距离是0，因此将满输入分配给这个容器。到下一个容器的距离是1，1-d=0，所以不分配给下一个容器，所有容器都是如此。假设某个方向值是45°，分配给第一个容器1/2，分配给第二个容器1/2。采用这种方法，为每个容器得到一个技术的比例分数，从而避免“边界”效应，即描述子方向的细微变化导致的为不同容器分配到同一个值的效应。
16个4×4区域得到16个直方图，将直方图的8个方向显示为一个小向量簇，其中每个向量的长度等于其对应容器的值。所以，一个描述子由4×4阵列组成，每个阵列包含8个方向值。在SIFT中，这描述子数据被组织为一个128维的向量。
3.归一化处理，将特征向量长度进行归一化处理，进一步去除光照的影响。
为了降低光照的影响，对特征向量进行了两阶段的归一化处理。首先，通过将每个分量除以向量范数，把向量归一化为单位长度。由每个像素值乘以一个常数所引起的图像对比度的变化，将以相同的常数乘以梯度，因此对比度的变化将被第一次归一化抵消。每个像素加上一个常量导致的亮度变化不会影响梯度值，因为梯度值是根据像素插值计算的。因此描述子对光照的仿射变换是不变的。然后，也可能出现摄像机饱和等导致的非线性光照变化。这类变化会导致某些梯度的相对幅值的较大变化，但它们几乎不会影响梯度方向。SIFT通过对归一化特征向量进行阈值处理，降低了较大梯度值的影响，使所有分量都小于实验确定的值0.2。阈值处理后，特征向量被重新归一化为单位向量。
#include
（1）对找到的极值点进行曲线插值拟合，并过滤adjustLocalExtrema()√
（2）计算方向直方图calcOrientationHist()√
（1）计算sift描述子calcSIFTDescriptor()
上面在讲原理的时候已经写了一些源码注释（打√），下面是剩下的注释。
2.2.1SIFT::createInitialImage()
其中σinit是第0层的尺度，σpre是被相机模糊后的尺度。在Lowe的论文使用了σ0=1.6，S=3。
为了尽可能多的保留原始图像信息，一般需要对原始图像进行扩大两倍采样，即升采样，从而生成一组采样图octave_1,此组采样图的第一层的模糊参数为：
也就是说，输入的图像已经是被模糊了（0.5），但你想要更模糊一点（1.6）的作为base。如果要先扩大图像尺寸的话，就先升采样后再用1.25的σ进行卷积。如果不升采样的话，就用sqrt(1.6*1.6-0.5*0.5)对原图做卷积，然后将结果作为base。
在前面的步骤中我们分别得到了关键点的位置（层，组，所在图像的横纵坐标），尺度，主方向信息。所以我们定位关键点到其所属的高斯金字塔中的那幅图中。用关键点周围(radius*2+1)*(radius*2+1)的邻域（橘色矩形框区域）来描述这个关键点，即计算邻域的幅值和幅角，其中
将这块区域分成d×d个小区域，SIFT中d取4，每个小区域的宽度是3σ，这里的σ是指相对于当前组第一层图片来说的。
2.校正旋转主方向，确保旋转不变性。
为了保证特征矢量的旋转不变性，要以特征点为中心，在附近邻域内将坐标轴旋转θ（特征点的主方向）角度，即将坐标轴旋转为特征点的主方向。如下图所示
旋转后邻域内像素的新坐标为：
3.建立三维直方图
在计算特征点描述符的时候，我们不需要精确知道邻域内所有像素的梯度幅值和幅角，我们只需要根据直方图知道其统计值即可。如图所示，三维直方图是由d*d*n个长宽为1的单位立方体组成的，高对应邻域像素幅角的大小，把360度分成8等分。立方体的底就是特征点的邻域区域，该区域被划分为4×4个子区域，邻域中的像素根据坐标位置，把它们归属到这16个子区域中的一个，再根据该像素幅角的大小，把它分到这8等分中的一份，这样每个像素点都能对应到其中的一个立方体内，三维直方图建立起来了。
正方体的中心代表着该正方体，但是落入正方体内的邻域像素不可能都在中心，因此我们需要对上面的梯度幅值做进一步的处理，根据它对中心点位置的贡献大小进行加权处理，即在正方体内，根据像素点相对于正方体中心的距离，对梯度幅值做加权处理。所以，三维直方图的值，也就是正方体的值需要下面四个步骤完成：
1）计算落入该正方体内的邻域像素的梯度幅值A
2）根据该像素相对于特征点的距离，对A进行高斯加权处理，得到B
3）根据该像素相对于它所在的正方体的中心的贡献大小，再对B进行加权处理，得到C
由于计算相对于正方体中心点的贡献大小略显繁琐,例如像素（0.25,0.25,0.25），中心点为（0.5,0.5,0.5），相对距离是（0.25,0.25,0.25），贡献权重为(1-0.25)*(1-0.25)*(1-0.25)。因此在实际应用中，我们需要经过坐标平移，把中心点平移到正方体的顶点上,这样只要计算正方体内的点对正方体的8个顶点的贡献即可。根据三线性插值法，对某个顶点的贡献值是以该顶点和正方体内的点为对角线的两个顶点，所构成的立方体的体积。也就是对8个顶点的贡献分别为：
经过上述处理后我们得到128维描述子，在OpenCV中，一共对描述子进行了两次归一化处理。第一次使用以下公式对描述子进行归一化处理，并对大于0.2的描述子进行截断处理：
假设经过第一次归一化后描述子用qi表示，第二次归一化操作是