sift
- 缘起
- 第一重：了解sift算子
- 1.1 基于局部特征的图像匹配
- 1.2 什么是sift算子
- 1.3 应用案例
- 1.3.1 参考链接
- 1.3.2 图像拼接与全景构建
- 1.3.2 图像配准
- 第二重：sift算子原理
- 2.1 图像的卷积滤波操作
- 2.2 多尺度空间
- 2.2.1 图像金字塔
- 2.2.2 高斯金字塔
- 2.2.3 高斯差分塔（ DOG算子）
- 2.3 特征提取
- 2.3.1 查找关键点(初步探查)
- 2.3.2 关键点定位
- 2.3.3 关键点方向分配
- 2.3.4 关键点特征描述(descriptors)
- step #1 确定计算描述子所需的图像区域
- step #2将坐标轴旋转为关键点的方向
- step #3 分配采样点
- step #4 计算每个种子点八个方向的梯度
- step #5 归一化处理
- step #6 描述子向量门限
- step #7 对特征描述向量排序
- 2.4 特征匹配查找
- 2.4.1 K-D tree
- 2.4.2 BBF
- 2.5 消除错配&变换矩阵
- 2.5.1 随机抽样一致算法（RANSAC)
- 2.5.2 RANSAC去除错配
- 2.6 改进的sift算子
- 2.6.1 PCA-SIFT
- 2.6.2 SURF
- 第三重：sift证明与源码分析
- 3.1 sift相关理论证明
- 3.1.1 消除边缘响应
- 3.1.2 Hessian矩阵
- 3.2 源码分析
- 3.3 sift缺点
- 3.1 sift相关理论证明
- 后记
缘起
学习sift算子(scale-invariant feature transform)的过程费时费力却受益匪浅。一方面，是因为sift算子本身包含内容丰富，另一方面，是因为学习过程中需要不断补充数值分析、数字图像处理等学科基础。虽然网上有很多朋友对sift算子写的不错，但是宏观上缺乏清晰的学习进阶思路，细节上对许多数学问题又描述不够。得益于Lowe的文章与学者对sift算子的讨论研究，我尝试总结一下对sift的理解。本文在兼顾通俗易懂的基础上，力求全面，希望成为一片完整概括和介绍sift算子的导论性文章。
第一重：了解sift算子
1.1 基于局部特征的图像匹配
基于局部特征的图像匹配（local feature-based image matching）算法在许多摄影测量、遥感和计算机视觉等领域有着重要应用，如图像配准（image registration）、变化检测(change detection)、图像融合（image fusion）和3D重构与建模(3D reconstruction and modeling)。局部特征可以定义为具有不同于其邻近区域的特性的图像模式或不同的结构。局部特征最重要的特性是鲁棒性，鲁棒性可以定义为特征对各种图像几何和光度变换的稳定性。除了特征鲁棒性之外，大多数局部特征检测都提取了以下三种类型的特征：
旋转不变特征(Rotation-invariant features)：旋转不变特征，也称为角点（corners and points），对应于二维图像中曲率较大或灰度变化较大的点。常用的角点特征检测算法有Harris（Harris and Stephens, 1988）、Hessian(Beaudet, 1978)、phase congruency(Kovesi, 2003)、FAST(Rosten et al., 2010)等。角点特征对平移和旋转不变，对光照变化和视点的小范围变化也具有一定的有鲁棒性。但是，它们对尺度变化敏感。
尺度不变特征(Scale-invariant features)：为了处理尺度变化，引入了尺度不变方法，利用尺度空间理论自动确定局部特征的位置和尺度(Lindeberg, 2008)。典型的有圆形特征。大多数尺度不变特征检测算子都是基于高斯二阶导数滤波器(Mainali et al., 2013)。著名的尺度不变特征检测算子有：automatic scale selection based on the scale-normalized Laplacian of Gaussian(Lindeberg,1998)，scale-invariant feature transform(SIFT) (Lowe,2004)，Harris-Hessian Laplace(Mikolajczyk and Schmid, 2004)， speeded-up robust features (SURF) (Bay et al., 2008)，scale-invariant feature operator(SFOP)(Forstner et al .,2009)，和binary robust invariant scalable keypoints (BRISK) (Leutenegger et al., 2011)。
仿射不变量特征（Affine-invariant features）：单一尺度的特征提取不能有效处理视点（viewpoint）变化带来的影响。为了在视点变化的大尺度图片中实现可靠的特征匹配，引入了仿射不变量方法(Mikolajczyk等人，2005)。在最初，大多数的仿射不变算法的输出形状为椭圆，如Harris/Hessian-Affine (Mikolajczyk and Schmid, 2004)，the intensity extremal-based region detector(IBR)(Tuytelaars and Van Gool, 2004)，and the salient region detector (Kadir et al., 2004)。然而，为了简单起见，其他类型提取的特征也被一个椭圆所替代。典型的有，irregular maximally stable extremal regions(MSER) (Matas et al.， 2004)和the parallelogram-shaped features of the edge-based region (EBR) detector (Tuytelaars and Van Gool, 2004)的平行四边形特征。
由于尺度不变特征和仿射不变特征提取通常需要进行特征尺度选择和仿射形状估计，所以其计算复杂度通常比点特征检测器算法高。局部特征通常与描述符/描述子（descriptors）相关联来描述和匹配它们。特征描述符通常是编码局部特征图像邻域属性的数值向量。有很多文献提出不同的描述符编码方法，其中SIFT (Lowe, 2004)、DAISY (Tola et al., 2010)、 LIOP (Wang et al., 2011)和crisp (Leutenegger et al., 2011)是局部特征描述符的几个值得注意的例子。通过对两幅图像的局部特征检测（生成特征点）和描述（生成描述子），利用欧氏距离或马氏距离等特定的相似度测度建立对应过程。
1.2 什么是sift算子
尺度不变特征转换(Scale-invariant feature transform或SIFT)是一种电脑视觉的算法用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量，此算法由 David Lowe在1999年所发表，2004年完善总结。其应用范围包含物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对。
此算法有其专利，专利拥有者为英属哥伦比亚大学。
局部影像特征的描述与侦测可以帮助辨识物体，SIFT 特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关。对于光线、噪声、些微视角改变的容忍度也相当高。基于这些特性，它们是高度显著而且相对容易撷取，在母数庞大的特征数据库中，很容易辨识物体而且鲜有误认。使用 SIFT特征描述对于部分物体遮蔽的侦测率也相当高，甚至只需要3个以上的SIFT物体特征就足以计算出位置与方位。在现今的电脑硬件速度下和小型的特征数据库条件下，辨识速度可接近即时运算。SIFT特征的信息量大，适合在海量数据库中快速准确匹配。
SIFT算法的特点有：
- SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性；
- 独特性（Distinctiveness）好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配；
- 多量性，即使少数的几个物体也可以产生大量的SIFT特征向量；
- 高速性，经优化的SIFT匹配算法甚至可以达到实时的要求；
- 可扩展性，可以很方便的与其他形式的特征向量进行联合。
SIFT算法可以解决的问题：目标的自身状态、场景所处的环境和成像器材的成像特性等因素影响图像配准/目标识别跟踪的性能。而SIFT算法在一定程度上可解决：
- 目标的旋转、缩放、平移（RST）
- 图像仿射/投影变换（视点viewpoint）
- 光照影响（illumination）
- 目标遮挡（occlusion）
- 杂物场景（clutter）
- 噪声
SIFT算法的实质是在不同的尺度空间上查找关键点(特征点)，并计算出关键点的方向。SIFT所查找到的关键点是一些十分突出，不会因光照，仿射变换和噪音等因素而变化的点，如角点、边缘点、暗区的亮点及亮区的暗点等。
1.3 应用案例
1.3.1 参考链接
C++-OpenCV封装
https://docs.opencv.org/3.1.0/d5/d3c/classcv_1_1xfeatures2d_1_1SIFT.html
Python-OpenCV封装；
https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_feature2d/py_sift_intro/py_sift_intro.html
Matlab里面有对surf算子的封装（我只知道2017b支持，2014a不支持。不太清楚具体从哪一版本加进去的）：
https://www.mathworks.com/help/gpucoder/examples/feature-extraction-using-surf.html
用户使用MATLAB对sift算子的封装：
https://www.mathworks.com/matlabcentral/fileexchange/29004-feature-points-in-image-keypoint-extraction
https://www.mathworks.com/matlabcentral/fileexchange/50319-sift-feature-extreaction
1.3.2 图像拼接与全景构建
为了构建图像全景图，我们将使用计算机视觉和图像处理技术，如:关键点检测和局部不变描述符;关键点匹配;RANSAC;和视角扭曲。
我们的全景拼接算法包括四个步骤:
step #1：检测关键点(DoG, Harris等)并从两个输入图像中提取局部不变描述符(SIFT, SURF等)。
step #2：匹配两个图像之间的描述符。
step #3：使用RANSAC算法，用我们的匹配特征向量估计一个单应性矩阵。
step #4：使用从步骤3得到的单应性矩阵应用扭曲变换。
我们在panorama.py文件中封装这四个步骤，在这个文件中，我们定义一个Stitcher的类来建立全景图。
Stitcher类依赖imutils包，可以使用下面语句在系统中安装imutils
pip install imutils下面介绍panorama.py
# import the necessary packages
import numpy as np
import imutils
import cv2
class Stitcher:
def __init__(self):
# determine if we are using OpenCV v3.X
self.isv3 = imutils.is_cv3()上面2-4行导入必要的包。在第6行定义了Stitcher类，构造函数通过调用imutils.is_cv3()方法检查OpenCV版本是否是v3.X。下一步，在Stitcher类中定义stitch()方法
def stitch(self, images, ratio=0.75, reprojThresh=4.0,
showMatches=False):
# unpack the images, then detect keypoints and extract
# local invariant descriptors from them
(imageB, imageA) = images
(kpsA, featuresA) = self.detectAnddescribe(imageA)
(kpsB, featuresB) = self.detectAndDescribe(imageB)
# match features between the two images
M = self.matchKeypoints(kpsA, kpsB,
featuresA, featuresB, ratio, reprojThresh)
# if the match is None, then there aren't enough matched
# keypoints to create a panorama
if M is None:
return None
# otherwise, APPly a perspective warp to stitch the images
# together
(matches, H, status) = M
result = cv2.warpPerspective(imageA, H,
(imageA.shape[1] + imageB.shape[1], imageA.shape[0]))
result[0:imageB.shape[0], 0:imageB.shape[1]] = imageB
# check to see if the keypoint matches should be visualized
if showMatches:
vis = self.drawMatches(imageA, imageB, kpsA, kpsB, matches,
status)
# return a tuple of the stitched image and the
# visualization
return (result, vis)
# return the stitched image
return resultstitch()方法只需要一个参数：images。images是由一个两张需要拼接图片构成的列表形式(两张图片必须按照从左到右的顺序)。其他参数有：ratio-用于David Lowe的比值测试，用于匹配特征；reprojThresh-RANSAC算法使用的最大摆动空间；showMatches-用于指示关键点匹配是否应该可视化的布尔值。
下面是stitch()方法调用的子方法。注意，这还是Stitcher类的方法。detectAndDescribe方法接受图像，然后检测关键点并提取局部不变描述符。在我们的实现中，我们使用了高斯差分塔(DoG)关键点探测器和SIFT特征提取器。
def detectAndDescribe(self, image):
# convert the image to grayscale
gray = cv2.cvtcolor(image, cv2.COLOR_BGR2GRAY)
# check to see if we are using OpenCV 3.X
if self.isv3:
# detect and extract features from the image
descriptor = cv2.xfeatures2d.SIFT_create()
(kps, features) = descriptor.detectAndCompute(image, None)
# otherwise, we are using OpenCV 2.4.X
else:
# detect keypoints in the image
detector = cv2.FeatureDetector_create("SIFT")
kps = detector.detect(gray)
# extract features from the image
extractor = cv2.DescriptorExtractor_create("SIFT")
(kps, features) = extractor.compute(gray, kps)
# convert the keypoints from KeyPoint objects to NumPy
# arrays
kps = np.float32([kp.pt for kp in kps])
# return a tuple of keypoints and features
return (kps, features)下面介绍matchKeypoints方法。用来匹配关键点，并生成单应性矩阵。
def matchKeypoints(self, kpsA, kpsB, featuresA, featuresB,
ratio, reprojThresh):
# compute the raw matches and initialize the list of actual
# matches
matcher = cv2.DescriptorMatcher_create("BruteForce")
rawMatches = matcher.knnMatch(featuresA, featuresB, 2)
matches = []
# loop over the raw matches
for m in rawMatches:
# ensure the distance is within a certain ratio of each
# other (i.e. Lowe's ratio test)
if len(m) == 2 and m[0].distance < m[1].distance * ratio:
matches.append((m[0].trainIdx, m[0].queryIdx))
# computing a homography requires at least 4 matches
if len(matches) > 4:
# construct the two sets of points
ptsA = np.float32([kpsA[i] for (_, i) in matches])
ptsB = np.float32([kpsB[i] for (i, _) in matches])
# compute the homography between the two sets of points
(H, status) = cv2.findHomography(ptsA, ptsB, cv2.RANSAC,
reprojThresh)
# return the matches along with the homograpy matrix
# and status of each matched point
return (matches, H, status)
# otherwise, no homograpy could be computed
return NonematchKeypoints函数需要四个参数：与第一张图像相关联的关键点和特征向量，然后是与第二张图像相关联的关键点和特征向量。还有David Lowe的比值检验变量和RANSAC重投影阈值。
将特性匹配在一起实际上是一个相对简单的过程。我们从两个图像中遍历描述符，计算距离，并找到每对描述符的最小距离。由于这是计算机视觉中的一种非常常见的实践，OpenCV的cv2.DescriptorMatcher_create内置函数，用来构造特征匹配器。BruteForce值表明计算机将从两个图像中详尽地计算所有特征向量之间的欧氏距离，并找到距离最小的描述符对。
使用knnMatch的调用使用k=2在两个特征向量集之间执行k- nn匹配(k=2表示返回每个特征向量的前两个匹配)。我们想要前两个匹配对而不仅仅是前一对的原因是我们需要应用ratio测试来进行假阳性匹配修剪。
计算了每一对描述符的rawMatches——但是有可能其中一些是假阳性，也就是说图像补丁实际上没有真的匹配。为了修剪这些假阳性匹配，我们单独遍历每个rawMatches并应用ratio测试，该测试用于确定高质量的特性匹配。比率一般设置在[0.7,0.8]范围内。（简单来说，就是对于image1的某个关键点，它跟image2的所有关键点去匹配。匹配的最好结果也就是最小distance0，必须比匹配的次最好匹配结果distance1*ratio还要小，才认为这两个关键点匹配得上）。
通过ratio测试获得匹配对之后，我们可以计算这两组关键点之间的单应性矩阵。在两组点之间计算单应性至少需要初始的四组匹配。为了获得更可靠的单应性估计，我们应该有比4个匹配点多得多的匹配点。
最后，介绍Stitcher类中的drawMatches方法用来可视化两幅图片的关键点对。
def drawMatches(self, imageA, imageB, kpsA, kpsB, matches, status):
# initialize the output visualization image
(hA, wA) = imageA.shape[:2]
(hB, wB) = imageB.shape[:2]
vis = np.zeros((max(hA, hB), wA + wB, 3), dtype="uint8")
vis[0:hA, 0:wA] = imageA
vis[0:hB, wA:] = imageB
# loop over the matches
for ((trainIdx, queryIdx), s) in zip(matches, status):
# only process the match if the keypoint was successfully
# matched
if s == 1:
# draw the match
ptA = (int(kpsA[queryIdx][0]), int(kpsA[queryIdx][1]))
ptB = (int(kpsB[trainIdx][0]) + wA, int(kpsB[trainIdx][1]))
cv2.line(vis, ptA, ptB, (0, 255, 0), 1)
# return the visualization
return vis该方法要求我们通过两个原始图像、两张图像对应的两组关键点，应用ratio检验后的关键点对，最后通过单应性计算提供的状态列表。使用这些变量，我们可以通过从第一张图片中的关键点N到第二张图片中的关键点M绘制一条直线来可视化内点。
最后还有一个stitch.py驱动脚本，用来简单地处理加载图像，调整它们的大小(以适配屏幕)，并构建全景图。
# import the necessary packages
from stitch.panorama import Stitcher
import imutils
import cv2
# load the two images and resize them to have a width of 400 pixels
# (for faster processing)
imageA = cv2.imread("untitled1.jpg")
imageB = cv2.imread("untitled2.jpg")
imageA = imutils.resize(imageA, width=400)
imageB = imutils.resize(imageB, width=400)
# stitch the images together to create a panorama
stitcher = Stitcher()
(result, vis) = stitcher.stitch([imageA, imageB], showMatches=True)
# show the images
cv2.imshow("Image A", imageA)
cv2.imshow("Image B", imageB)
cv2.imshow("Keypoint Matches", vis)
cv2.imshow("Result", result)
cv2.waitKey(0)提示：
1. 如果OpenCV出现 ‘module’ object has no attribute ‘xfeatures2d’问题，点击这里；
2. 调用函数时，两幅图片从左必须服从从左到右的顺序。
1.3.2 图像配准
import cv2
print("sift is so interesting")至此，介绍完sift算子第一重。对于只关心sift算子功能和用法的朋友已经足够，不必再更近第二重，深入了解sift算子的原理了。
第二重：sift算子原理
2.1 图像的卷积滤波操作
2.1.1滤波、模板与卷积核
参考冈萨雷斯数字图像处理滤波等章节内容。
2.1.2高斯模糊
SIFT算法是在不同的尺度空间上查找关键点，而尺度空间的获取需要使用高斯模糊来实现，Lindeberg等人已证明高斯卷积核是实现尺度变换的唯一变换核，并且是唯一的线性核。本节先介绍高斯模糊算法。
高斯模糊是一种图像滤波器，它使用正态分布(高斯函数)计算模糊模板，并使用该模板与原图像做卷积运算，达到模糊图像的目的。N维空间正态分布方程为：
其中，
在二维空间中，这个公式生成的曲面的等高线是从中心开始呈正态分布的同心圆，如下图所示。分布不为零的像素组成的卷积矩阵与原始图像做变换。每个像素的值都是周围相邻像素值的加权平均。原始像素的值有最大的高斯分布值，所以有最大的权重，相邻像素随着距离原始像素越来越远，其权重也越来越小。这样进行模糊处理比其它的均衡模糊滤波器更高地保留了边缘效果。
理论上来讲，图像中每点的分布都不为零，这也就是说每个像素的计算都需要包含整幅图像。在实际应用中，在计算高斯函数的离散近似时，在大概
根据
| 0.0165 | 0.0297 | 0.0361 | 0.0297 | 0.0165 |
| 0.0297 | 0.0534 | 0.0649 | 0.0534 | 0.0297 |
| 0.0361 | 0.0649 | 0.0789 | 0.0649 | 0.0361 |
| 0.0297 | 0.0534 | 0.0649 | 0.0534 | 0.0297 |
| 0.0165 | 0.0297 | 0.0361 | 0.0297 | 0.0165 |
使用高斯模糊对图像处理结果如下图所示。
下面的MATLAB代码用来演示高斯滤波的使用方法：
代码运行环境：thinkpad e480；windows10 教育版；MATLAB2017b
clear all;
clc;
%----------------------------------------------
%对图像进行高斯滤波，并显示图像
%----------------------------------------------
%读进图像
[filename, pathname] = uigetfile({'*.jpg'; '*.bmp'; '*.gif'; '*.png' }, '选择图片');
%没有图像
if filename == 0
return;
end
Image = imread([pathname, filename]);
[m, n, z] = size(Image);
%转换为灰度图
if z>1
Image = rgb2gray(Image);
end
sigma = 1;
gausFilter = fspecial('gaussian', [5,5], sigma);
gaus= imfilter(Image, gausFilter, 'replicate');
%显示图像-----------------------
figure(1)
subplot(1,2,1);
imshow(Image);
title('原图像');
subplot(1,2,2);
imshow(gaus);
title('滤波后');2.2 多尺度空间
2.2.1 图像金字塔
图像金字塔是一种以多分辨率来解释图像的结构，通过对原始图像进行多尺度像素采样的方式，生成N个不同分辨率的图像。把具有最高级别分辨率的图像放在底部，以金字塔形状排列，往上是一系列像素（尺寸）逐渐降低的图像，一直到金字塔的顶部只包含一个像素点的图像，这就构成了传统意义上的图像金字塔。
获得图像金字塔一般包括两个步骤：
1. 利用低通滤波器平滑图像
2. 对平滑图像进行抽样（采样）有两种采样方式——上采样（分辨率逐级升高）和下采样（分辨率逐级降低）
FLAG2： 以后补充上采样、下采样内容。
2.2.2 高斯金字塔
高斯金字塔是在sift算子中提出来的概念，首先高斯金字塔并不是一个图像金字塔，而是有很多组（Octave）金字塔构成，并且每组金字塔都包含若干层（Interval）。
高斯金字塔构建过程：
对于参数σ，在Sift算子中取的是固定值1.6。
- 先将原图像扩大一倍之后作为高斯金字塔的第1组第1层，将第1组第1层图像经高斯卷积（其实就是高斯平滑或称高斯滤波）之后作为第1组金字塔的第2层。
将
σ k " role="presentation">,等到一个新的平滑因子σ = k ∗ σ 如此这般重复，最后得到L层图像，在同一组中，每一层图像的尺寸都是一样的，只是平滑系数不一样。它们对应的平滑系数分别为：
0 " role="presentation">,σ " role="presentation">，k σ k 2 σ k 3 σ k ( L − 2 ) σ " role="presentation">。将第1组倒数第三层（也就是从上往下数第三层）图像作比例因子为2的降采样，得到的图像作为第2组的第1层，然后对第2组的第1层图像做平滑因子为σ的高斯平滑，得到第2组的第2层，就像步骤2中一样，如此得到第2组的L层图像，同组内它们的尺寸是一样的，对应的平滑系数分别为：
0 ， σ ， k σ ， k 2 σ , k 3 σ … … k ( L − 2 ) σ
这样反复执行，就可以得到一共
2.2.3 高斯差分塔（ DOG算子）
在实际计算时，使用高斯金字塔每组中相邻上下两层图像相减，得到高斯差分塔（Difference of Gaussian， DOG算子），如下图所示，进行极值检测。
表示为：
符号表示：
在上述尺度空间中，
其中
在最开始建立高斯金字塔时，要预先模糊输入图像来作为第0个组的第0层的图像，这时相当于丢弃了最高的空域的采样率。因此通常的做法是先将图像的尺度扩大一倍来生成第-1组。我们假定初始的输入图像为了抗击混淆现象，已经对其进行
取式(2-3)中的k为组内总层数的倒数，即
在构建高斯金字塔时，组内每层的尺度坐标按如下公式计算：
其中
该组内尺度在方向分配和特征描述时确定采样窗口的大小。
由上，式(2-3)可记为
下图为构建DOG金字塔的示意图，原图采用128*128的jobs图像，扩大一倍后构建金字塔。
2.3 特征提取
2.3.1 查找关键点(初步探查)
为了寻找尺度空间的极值点，每一个采样点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。如下图所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。 一个点如果在DOG尺度空间本层以及上下两层的26个领域中是最大或最小值时，就认为该点是图像在该尺度下的一个特征点,如图所示。
在极值比较的过程中，每一组图像的首末两层是无法进行极值比较的，为了满足尺度变化的连续性，我们在每一组图像的顶层继续用高斯模糊生成了 3 幅图像，高斯金字塔有每组S+3层图像。DOG金字塔每组有S+2层图像.（如下图所示）
为了满足尺度变化的连续性：假设s=3，也就是每个塔里有3层，则k=21/s=21/3，那么按照上图可得Gauss Space和DoG space 分别有3个（s个）和2个（s-1个）分量，在DoG space中，1st-octave两项分别是σ,kσ; 2nd-octave两项分别是2σ,2kσ;由于无法比较极值，我们必须在高斯空间继续添加高斯模糊项，使得形成σ,kσ,k2σ,k3σ,k4σ这样就可以选择DoG space中的中间三项kσ,k2σ,k3σ（只有左右都有才能有极值），那么下一octave中（由上一层降采样获得）所得三项即为2kσ,2k2σ,2k3σ，其首项2kσ=24/3。刚好与上一octave末项k3σ=23/3尺度变化连续起来，所以每次要在Gaussian space添加3项，每组（塔）共S+3层图像，相应的DoG金字塔有S+2层图像。
当然这样产生的极值点并不全都是稳定的特征点，因为某些极值点响应较弱，而且DOG算子会产生较强的边缘响应。
注意：使用Laplacian of Gaussian能够很好地找到找到图像中的兴趣点，但是需要大量的计算量，所以使用Difference of Gaussian图像的极大极小值近似寻找特征点。DOG算子计算简单，是尺度归一化的LoG算子的近似。
2.3.2 关键点定位
以上方法检测到的极值点是离散空间的极值点，以下通过拟合三维二次函数来精确确定关键点的位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点(因为DoG算子会产生较强的边缘响应)，以增强匹配稳定性、提高抗噪声能力。
离散空间的极值点并不是真正的极值点，下图显示了二维函数离散空间得到的极值点与连续空间极值点的差别。利用已知的离散空间点插值得到的连续空间极值点的方法叫做子像素插值（Sub-pixel Interpolation）。
为了提高关键点的稳定性，需要对尺度空间DoG函数进行曲线拟合。利用DoG函数在尺度空间的Taylor展开式(拟合函数)为：
其中，
对应极值点，方程的值为：
其中，
2.3.3 关键点方向分配
为了使描述符具有旋转不变性，需要利用图像的局部特征(关键点邻域像素的梯度方向分布特性)给每一个关键点分配一个基准方向。对于在DOG金字塔中检测出的关键点，采集其所在高斯金字塔图像
在完成关键点的梯度计算后，使用直方图统计邻域内像素的梯度和方向。梯度直方图的范围是0～360度，其中每10度一个柱，总共36个柱；或者每45度一个柱，总共8个柱。如下图所示，直方图的峰值则代表了该关键点处邻域梯度的主方向，即作为该关键点的方向。下图所示为八个方向的梯度直方图。
邻域内的像素距中心点越远，其对直方图的贡献也响应也越小，Lowe论文中提到使用高斯函数对直方图进行平滑。按Lowe的建议，梯度的模值
梯度直方图的峰值代表了该特征点处邻域梯度的方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向作为该关键点的辅方向。因此，对于同一梯度值的多个峰值的关键点位置，在相同位置和尺度将会有多个关键点被创建但方向不同。仅有15％的关键点被赋予多个方向，但可以明显的提高关键点匹配的稳定性。实际编程实现中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点，并且，离散的梯度方向直方图要进行插值拟合处理，来求得更精确的方向角度值，检测结果下图所示。
至此，将检测出的含有位置、尺度和方向的关键点即是该图像的SIFT特征点。
2.3.4 关键点特征描述(descriptors)
通过以上步骤，对于每一个关键点，拥有三个信息：位置、尺度以及方向。接下来就是为每个关键点建立一个描述符，用一组向量将这个关键点描述出来，使其不随各种变化而改变，如光照变化、视角变化等。这个描述子不但包括关键点，也包含关键点周围对其有贡献的像素点，并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。
SIFT描述子是关键点邻域高斯图像梯度统计结果的一种表示。通过对关键点周围图像区域分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性。
Lowe建议描述子使用在关键点尺度空间内
step #1 确定计算描述子所需的图像区域
特征描述子与特征点所在的尺度有关，因此，对梯度的求取应在特征点对应的高斯图像上进行。将关键点附近的邻域划分为d*d(Lowe建议d=4)个子区域，每个子区域做为一个种子点，每个种子点有8个方向。每个子区域的大小与关键点方向分配时相同，即每个区域有个
计算结果四舍五入为整数。
step #2将坐标轴旋转为关键点的方向
将坐标轴旋转为关键点的方向，以确保旋转不变性，如下图所示。
旋转后邻域内采样点的新坐标为：
step #3 分配采样点
将邻域内的采样点分配到对应的子区域内，将子区域内的梯度值分配到8个方向上，计算其权值。
旋转后的采样点坐标在半径为
旋转后的采样点(x’,y’)落在子区域的下标为
Lowe建议子区域的像素的梯度大小按
其中
step #4 计算每个种子点八个方向的梯度
如图6.3所示，将由式(6-3)所得采样点在子区域中的下标(图中蓝色窗口内红色点)线性插值，计算其对每个种子点的贡献。如图中的红色点，落在第
其中
step #5 归一化处理
如上统计的
step #6 描述子向量门限
非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此设置门限值(向量归一化后，一般取
step #7 对特征描述向量排序
按特征点的尺度对特征描述向量进行排序。
到这里sift算子运算流程已经基本交代完毕，但是在实际应用过程中又常以sift算子为基础，将生成的特征进行匹配并生成变换矩阵，这里一并介绍。
2.4 特征匹配查找
SIFT特征匹配算法，和数据库查、图像检索本质上是同一个问题，都可以归结为一个通过距离函数在高维矢量空间进行相似性检索的问题，如何快速而准确地找到查询点的近邻，不少人提出了很多高维空间索引结构和近似查询的算法。
一般说来，索引结构中相似性查询有两种基本的方式：
- 一种是范围查询，范围查询时给定查询点和查询距离阈值，从数据集中查找所有与查询点距离小于阈值的数据；
- 另一种是K近邻查询，就是给定查询点及正整数K，从数据集中找到距离查询点最近的K个数据，当K=1时，它就是最近邻查询。
同样，针对特征点匹配也有两种方法：
- 最容易的办法就是线性扫描，也就是我们常说的穷举搜索，依次计算样本集E中每个样本到输入实例点的距离，然后抽取出计算出来的最小距离的点即为最近邻点。此种办法简单直白，但当样本集或训练集很大时，它的缺点就立马暴露出来了。举个例子，在物体识别的问题中，可能有数千个甚至数万个SIFT特征点，而去一一计算这成千上万的特征点与输入实例点的距离，明显是不足取的。
- 另外一种，就是构建数据索引，因为实际数据一般都会呈现簇状的聚类形态，因此我们想到建立数据索引，然后再进行快速匹配。索引树是一种树结构索引方法，其基本思想是对搜索空间进行层次划分。根据划分的空间是否有混叠可以分为Clipping和overlapping两种。前者划分空间没有重叠，其代表就是K-D树；后者划分空间相互有交叠，其代表为R树。
在K-D tree上查找KNN使用BBF（多说一句，在R数上查找也有相应的算法，如BAB算法是最早的基于R数的静态KNN查询算法）。
2.4.1 K-D tree
参考文件第七章，讲的挺详细的，这里就不重复了。
2.4.2 BBF
https://blog.csdn.net/v_july_v/article/details/8203674
2.5 消除错配&变换矩阵
2.5.1 随机抽样一致算法（RANSAC)
随机抽样一致算法（RANdom SAmple Consensus，RANSAC）是一种采用迭代方式从一组包含离群的被观测数据中估算出数学模型参数的一种方法。
RANSAC算法包含两个基本假设：样本中包含正确数据(inliers，符合模型的数据)和异常数据(Outliers，不符合模型的数据)，即数据集中含有噪声。这些异常数据可能是由于错误的测量、错误的假设、错误的计算等产生的；同时RANSAC也假设， 给定一组正确的数据，存在可以计算出符合这些数据的模型参数的方法。
RANSAC算法基本思想描述 (仔细看，这里能看明白)：
- 考虑一个最小抽样集的势为n的模型(n为初始化模型参数所需的最小样本数)和一个样本集P，集合P的样本数#(P)>n，从P中随机抽取包含n个样本的P的子集S初始化模型M；
- 余集SC=P\S中与模型M的误差小于某一设定阈值t的样本集以及S构成S*。S*认为是内点集，它们构成S的一致集(Consensus Set)；
- 若#(S*)≥N，认为得到正确的模型参数，并利用集S*(内点inliers)采用最小二乘等方法重新计算新的模型M*；重新随机抽取新的S，重复以上过程。
- 在完成一定的抽样次数后，若未找到一致集则算法失败，否则选取抽样后得到的最大一致集判断内外点，算法结束。
RANSAC算法经常用于计算机视觉，例如同时求解相关问题与估计立体摄像机的基础矩阵，在图像拼接时求变换矩阵的时候。
想要了解更多RANSAC方面的细节，建议参考Marco Zuliani的Ransac for Dummies。此外，MATLAB官方教程也提供了RANSAC部分应用场景。
2.5.2 RANSAC去除错配
简单介绍完RANSAC基本原理，这里介绍RANSAC在SIFT算子里面的具体应用。RANSAC算法在SIFT特征筛选中的主要流程是：
- 从样本集中随机抽选一个RANSAC样本，即4个匹配点对
- 根据这4个匹配点对计算变换矩阵M（单应性矩阵）
- 根据样本集，变换矩阵M，和误差度量函数计算满足当前变换矩阵的一致集consensus，并返回一致集中元素个数
- 根据当前一致集中元素个数判断是否最优(最大)一致集，若是则更新当前最优一致集
- 更新当前错误概率p，若p大于允许的最小错误概率则重复1至4继续迭代，直到当前错误概率p小于最小错误概率
单应性矩阵：参考这里。资料有点长，看相关的知识就行了。
print("真的不能再展开了")
print("跟sift算子框架关系不大且具有严密推理过程的内容我会放上参考链接")2.6 改进的sift算子
2.6.1 PCA-SIFT
2.6.2 SURF
print("sift算子讲完了，根据所学知识举一反三...")
print("真不是我懒")此为sift算子第二重，第二重境界已经满足绝大部分人对sift算子运算过程的好奇心，然而在部分问题的理论证明上，并没有详细给出。这些问题将在第三重全部得到解答。但是进入第三重境界之前，必须提醒读者必须具有一定的数字图像处理基础与一定编程能力（这些内容可以在网络上获取到），不然你会跟我一样，陷入公式与程序的恐慌之中。
第三重：sift证明与源码分析
3.1 sift相关理论证明
3.1.1 消除边缘响应
一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。因此DOG算子会产生较强的边缘响应，需要剔除不稳定的边缘响应点。
首先，获取特征点处的Hessian矩阵，主曲率就是通过一个2x2 的Hessian矩阵H求出：
分别表示矩阵
Hessian矩阵是由采样点相邻差估计得到的，在下一节中说明。
D的主曲率和H的特征值成正比，令为α最大特征值，β为最小的特征值，则公式的值在两个特征值相等时最小，随着的增大而增大。值越大，说明两个特征值的比值越大，即在某一个方向的梯度值越大，而在另一个方向的梯度值越小，而边缘恰恰就是这种情况。所以为了剔除边缘响应点，需要让该比值小于一定的阈值，因此，为了检测主曲率是否在某域值r下，只需检测
上式成立时，关键点保留，反之剔除。
在Lowe的文章中，取r＝10。
3.1.2 Hessian矩阵
在数学中, 海森矩阵(Hessian matrix或Hessian)是一个自变量为向量的实值函数的二阶偏导数组成的方块矩阵, 此函数如下：
如果
其中，
(也有人把海森定义为以上矩阵的行列式)
Hessian矩阵有着广泛的应用，如在牛顿方法、求极值以及边缘检测、消除边缘响应等方面的应用。一个Hessian矩阵涉及到很多数学相关的知识点，如泰勒公式、极值判断、矩阵特征值及特征向量、二次型等。关于Hessian矩阵的应用在此不一一说明，这里只介绍利用有限差分法计算Hessian矩阵的方法。
利用泰勒公式证明有限差分法，参见维基百科
在图像处理中，取h=1，在上图中，将像素点位置0处的基本中点导数公式如下：
3.2 源码分析
print("source code analysis")3.3 sift缺点
SIFT在图像的不变特征提取方面拥有无与伦比的优势，却并非完美，其仍然存在：
- 实时性不高
- 有时特征点较少
- 对边缘光滑的目标无法准确提取特征点
等缺点。对模糊的图像和边缘平滑的图像，检测出的特征点过少，对圆更是无能为力。近来不断有人改进，著名的有Speeded Up Robust Features （SURF）和Colored SIFT （CSIFT） 等。
后记
sift算子中使用的多尺度空间，在图像处理中具有重要地位。其思想与小波分析理论的精髓多分辨率分析不谋而合，使我们能将sift算子与小波变换一定程度联系起来。对多尺度空间感兴趣的朋友，可以参考小波变换相关资料进行学习。
本文从2018年8月22日开始着手编辑，也是我第一次尝试写博客。本打算在博士入学前完工，后来偷工减料居然提前。这篇文章想尽可能全面地介绍sift算子，又想尽可能降低阅读门槛，所以在框架和内容上不断进行修改，希望对大家学习提高有所帮助。
最后，感谢诸多朋友的文章及著作，让我有机会在其基础上对sift算子进行总结、深入。由于本人能力有限，加上时间仓促，文章有任何不妥之处敬请广大读者随时批评、指正。感谢！e-mail： [email protected]
相关阅读
IP地址是以网络号和主机号来表示网络上的主机的，只有在一个网络号下的计算机之间才能“直接”互通，不同网络号的计算机要通过网关（Ga