SIFT特征详解

大名鼎鼎的SIFT特征对于每个做图像的人而言都不陌生，但是对于其细节却不是每个人都能讲清楚。基于最近研究Colmap中的稀疏重建的过程，准备把SIFT特征的细节盘一下，记录下来方便以后查阅。

SIFT的全称是Scale Invariant Feature Transform，它的几个特点是

尺度不变性
旋转不变性
不受亮度影响

这几个特点我们在下面的介绍中将具体解释。对于一个特征而言，主要有两部分，一部分是特征点，另一部分则是特征点的描述子。简单来说，就是说要找出图像中有个性的一些点，并且想办法独一无二地描述出来。

特征点的确定

SIFT特征点的获得需要构建高斯金字塔，它以输入图像为蓝本，构造不同分辨率的图像（octaves），针对某一个分辨率的图像，再利用不同尺度的高斯模糊获得一系列图像（sublevels）。在同一分辨率下，相邻的高斯模糊图像相减得到高斯差分图（Difference of Gaussian, DoG）。DoG图的极值点即对应了特征点。由于使用了不同分辨率的图像，不同尺度下的特征都可以找到。又由于使用了不同level的高斯模糊，则可以找到不同清晰程度下的特征点。最后，将不同分辨率不同高斯尺度找到的特征点映射回原图大小对应的图像坐标，即可得到所有的特征点。SIFT的尺度不变性就是由高斯金字塔带来的，因为它涵盖了各个尺度下的信息。

高斯模糊

对于一个图像 $I(x,y)$ 而言，其高斯模糊的图像可以表示为
$L(x,y,\sigma) = G(x,y,\sigma) \ast I(x,y)$
其中，
$G(x,y,\sigma)=\frac{1}{2\pi \sigma^2}\mathrm{e}^{-\frac{x^2+y^2}{2\sigma^2}}$
一般认为，高斯函数的值主要集中在 $\pm 3\sigma$ 内，因此在真正计算时，只需要采 $(6\sigma+1)\times(6\sigma+1)$ 的窗口。

金字塔的构建

图片 1.png

Octave数目的确定

输入图像为 $n\times n$ ，则对应Octave 0的图像大小为 $n\times n$ 。一般而言，会upsample到 $2n\times 2n$ ，形成Octave -1。这是因为Octave 0本身就需要做高斯模糊，这样输入图像就丢失了最清晰的状态，所以需要把它upsample再模糊，以期有一张图能保持原有的清晰。在Colmap的实现中（VLFeat或SiftGPU），设定了最小的Octave至少有 $16\times 16$ 个像素，于是Octave的数目便可以决定了
$n_o = \lfloor \log_2 n \rfloor - 2$

Sublevel 数目的确定

Sublevel数目是由参数octaveResolution 来决定的，假设octaveResolution为 $n_s$ ，则Sublevel的数目为 $n_s+3$ 。具体的原因我们之后会介绍。Sublevel的index也是从-1开始，到 $n_s+1$ 结束，总共 $n_s+3$ 个。Colmap中， $n_s = 3$ 。

每一个图像对应的模糊尺度 $\sigma$

首先会定义一个初始尺度 $\sigma_0$ ，该尺度对应的图像为Octave 0，Sublevel 0。其余图像的模糊尺度定义为
$\sigma(o,s) = \sigma_0 * 2^{o+\frac{s}{n_s}}$
在Colmap中，初始尺度为 $\sigma_0=1.6\times2^{\frac{1}{n_s}}$ 。容易看出，这样设置以后，Octave 0，Sublevel -1的图对应的尺度为 $\sigma = 1.6$ ，Octave -1， Sublevel -1的图像对应的尺度为 $\sigma = 0.8$ ，这与原始论文是一致的。根据上述公式，还可以发现 $\sigma(o,s) = \sigma(o-1, s+n_s)$ ，即Octave $o$ , Sublevel $s$ 的模糊尺度与Octave $o-1$ , Sublevel $s+n_s$ 的模糊尺度是一样的，而两张图差的只是一个resize。因此，在代码实现时，模糊只需要做一次，通过resize或者upsample操作就可以得到另一张图了，这样可以节省运算时间。

高斯差分图

对于每个固定的Octave，将相邻Sublevel的高斯模糊图相减即可得到高斯差分图。

DoG.png

请注意到，如果Sublevel的图像数量为 $n_s+3$ ，则DoG图像的数量为 $n_s+2$ 。DoG图可表示为
$D(x,y,k\sigma) = L(x,y,k\sigma) - L(x,y,(k-1)\sigma)$

极值的寻找

DoG图像空间的极值就是我们要寻找的极值点。对于DoG图像中的每个像素，只需要是它自身 $3\times3\times3$ 的邻域（本身图像 $3\times3$ 的邻域，上一个Sublevel图像 $3\times 3$ 的邻域，下一个Sublevel图像 $3\times 3$ 的邻域）的极大值，或者极小值，即是我们需要寻找的极值点。换言之，就是该像素的值比其他26个像素的值都大或者都小。

extreme.png

极值将不在最上面一层和最下面一层出现，因此极值只能出现在中间 $n_s$ 张图像上，这也就是octaveResolution的由来。

sublevel.png

极值的调整

极值的亚像素定位

由上述方法找出来的极值点显然都是整数的像素位置，为了更精确地得到极值的位置，我们需要利用Taylor展开确定极值的亚像素位置。我们以一维函数为例来说明。
$f(x+\delta x) = f(x) + f'(x)\cdot\delta x+\frac{1}{2}f''(x)\delta x^2$
对于极值位置 $x+\delta x$ ，有 $f'(x+\delta x) = 0$ ，上式两边对 $\delta x$ 求导，即有
$\delta x = - \frac{f'(x)}{f''(x)}$
代回即有
$f(x+\delta x) = f(x) + \frac{1}{2}f'(x)\delta x$
具体的实现中，上述过程迭代了5次，以找到更准确的极值位置。当 $\delta x$ 的绝对值大于0.5时，极值的位置将得到修改。

极值的筛选

对于DoG图像中绝对值较小的极值，本身容易受噪声干扰而不稳定，因此首先需要筛掉那些 $\|D(x,y,\sigma)\|$ 较小的极值点，Colmap中设定的阈值为0.01.

另外，应该尽量筛掉边缘的点，因为他们只有一个方向有极值，容易导致无法唯一地描述，影响后续匹配。边缘点可以通过Hessian矩阵的特征值的比值来找到，即边缘点Hessian矩阵的两个特征值会相差很大。对于Hessian矩阵
$H= \left(\begin{matrix} D_{xx} & D_{xy}\\ D_{xy} & D_{yy} \end{matrix}\right)$
直接计算其特征值并不明智，假设其特征值大的为 $\alpha$ ，小的为 $\beta$ ，我们利用
$\mathrm{Tr}(H) = D_{xx} + D_{yy} = \alpha + \beta, \det(H) = D_{xx}*D_{yy} - D_{xy}^2 = \alpha\beta$
可以得到
$\frac{\mathrm{Tr}^2(H)}{\det(H)} = \frac{(\alpha+\beta)^2}{\alpha\beta}= \frac{\left(\frac{\alpha}{\beta} + 1\right)^2}{\frac{\alpha}{\beta}}$
其中， $\mathrm{Tr}$ 表示矩阵的迹，而 $\det$ 表示矩阵的行列式。于是，我们得到了一个简单的计算 $\alpha$ 与 $\beta$ 比值的方法。对 $\frac{\alpha}{\beta}$ 设置一个阈值，即可筛掉一些边缘点。原论文中设置的阈值是10。在SiftGPU中，这个阈值也是10，但是它是针对 $\frac{\mathrm{Tr}^2(H)}{\det(H)}$ 来设置的。而在VLFeat中，这个阈值10又是针对公式
$\left(0.5*\frac{\mathrm{Tr}^2(H)}{\det(H)}-1\right) + \sqrt{\max\left(0.25*\frac{\mathrm{Tr}^2(H)}{\det(H)}, 0\right)*\frac{\mathrm{Tr}^2(H)}{\det(H)}}$
来设置的。如果按 $\frac{\alpha}{\beta}=10$ 来计算，SiftGPU中的阈值应该设成12.1，而VLFeat中的阈值应该设成11.1左右，总之这里各个实现之间是有些差异的。

特征点的尺度

经过极值的调整以后，每个特征点将会有一个亚像素的位置 $(x,y,z)$ ，其中 $x,y$ 对应其所在octave的像素坐标系，而 $z$ 对应某个中间的Sublevel图像index。我们可以重新计算该特征点的尺度为
$\sigma(o,z) = \sigma_0 2^{o+\frac{z-1}{n_s}}$

确定特征点方向

对于剩下的特征点，下一个操作是计算一个主方向，得到主方向以后，所有的特征点都可以旋转到主方向为 $+x$ 轴，以计算描述子。这也是Sift特征旋转不变性的原因。对于每个特征点，主方向对应的是邻域范围内梯度响应最大的那个方向。

梯度图的选取

首先我们需要确定梯度图是根据哪个尺度 $(o,s)$ 下的DoG图计算的，这个是由特征点的尺度来决定的，也就是说特征点是在哪个尺度下得到的，就在哪个尺度上计算梯度图。

邻域大小的确定

邻域为一个圆，其半径大小确定为 $3\sigma$ ，其中 $\sigma$ 对应该特征点的尺度。

角度方向直方图的生成

为了确定梯度响应最大的那个方向，我们需要生成一个梯度的角度方向直方图。我们将角度分为36份，每一份对应10度。在邻域范围内，计算每个像素的梯度幅值及角度，按照角度将该像素放置于哪个区间。在计算直方图时，我们累加的不是像素的个数，而是该像素的幅值，同时还会给其一个权重。这个权重值是以特征点为中心的一个高斯分布，其标准差为 $1.5\sigma$ ，其中 $\sigma$ 为该特征点的尺度。按照这个操作，我们最终可以得到一个角度方向直方图。

在得到初始的直方图之后，实现中还对其进行了6次平滑，主要是每个histogram值用其左右及自身的平均来代替。

histogram.png

主方向的确立

根据直方图，我们要得到最终的角度。具体的做法是，先确定直方图最高的bin，然后再利用左右的bin的幅值进行插值。比如，直方图最高的bin值为 $b$ ，对应的角度下限为 $\alpha$ ，其左右分别为 $a,c$ ，那么最终的角度为
$\alpha + 0.5*\frac{c-a}{2b-a-c} + 0.5$
另外，为了提高鲁棒性，我们不只使用了最高bin，而是把最高bin0.8倍以内的高值都记录下来了，当做一个方向。实现上，VLFeat中限制最多4个方向，而SiftGPU中只提供最多两个方向。多余的方向在具体的实现时，可以把特征点复制一份，赋给其新的方向即可。原始论文中提到，只有15%左右的特征点会拥有多个方向，但却能有效提高匹配的鲁棒性。

在SiftGPU中，为了减少显存空间并加速，它使用了一个float4个字节高位和低位分别存储两个方向，即每个方向对应2个字节，这时角度的精度为360/65535度。

特征点的输出

每一个特征点最终的输出包括 $(x,y,\sigma,\theta)$ ，其中 $x,y$ 位置将被scale到输入图像大小的像素空间，而 $\sigma$ 代表了该特征点的尺度， $\theta$ 对应该特征点的主方向。另外，还可以针对特征点进行排序，依据是Octave和Sublevel的index，从小到大排序。

特征点描述子

描述子的生成与计算角度方向时有些类似，思想是在特征点对应的尺度图像上，先将图像按主方向旋转，即将主方向转至 $+x$ 轴方向。再将旋转后的邻域范围等划分成一个 $4\times4$ 的子区域，在每个子区域上，我们计算其角度方向直方图，这时我们只采用8个bin。于是，对于每个特征点，我们将得到一个 $4\times4\times8=128$ 维的描述向量。

rotation.jpg

邻域大小的确定

实际实现中，并不需要真的对图像进行旋转，只需要把对应像素的坐标和梯度角度做相应变化即可。因此，在计算中，仍可使用原图，这里讨论的邻域大小也是相对原图而言。

首先我们期望每个子区域的大小为 $3\sigma$ ，其中 $\sigma$ 为特征点的尺度，那么至少需要邻域的边长为 $3\sigma \times 4$ 。为了插值的方便，一般我们会将其扩大一点，设置成 $3\sigma\times 5$ 。参考上图的旋转模式，要在旋转后保证邻域边长有 $3\sigma\times 5$ ，则需要在原图至少边长为 $3\sigma\times 5 \times \sqrt{2}$ ，于是邻域的半径就可以确定为

$\lfloor \frac{15 \sqrt{2}\sigma}{2}+0.5\rfloor$
在邻域范围内的每个点，我们将其旋转后，如果能落在特征点 $\pm 6\sigma$ 以内，则它可以对角度直方图有贡献。

角度方向直方图的计算

对于邻域中的每个点，我们先计算其在以特征点为中心的像素坐标系下的坐标 $(x,y)$ ，再将其旋转到主方向，得到主方向下的坐标 $(x',y')$ 。根据旋转公式有
$\left(\begin{matrix} x'\\ y' \end{matrix}\right) = \left(\begin{matrix} \cos\theta & -\sin\theta\\ \sin\theta & \cos\theta \end{matrix}\right) \left(\begin{matrix} x\\ y \end{matrix}\right)$
从 $(x',y')$ 的位置，我们通过计算 $\left(\lfloor\frac{x'}{3\sigma}\rfloor, \lfloor\frac{y'}{3\sigma}\rfloor\right)$ 可以推断出它属于哪个子区域。

division.png

实现中，它的梯度角度不只对它归属的子区域有贡献，还对旁边的一些子区域有贡献，具体哪些子区域有贡献，可以参考下图。图中的绿色区域代表 $4\times 4$ 的子区域，而蓝色区域代表更大的邻域区域 $(5\times5)$ ，假设一个点 $(x',y')$ 落在带阴影的蓝色区域内，该区域与四个子区域相交，分别是 $(-2,1), (-1,1), (-2,0), (-1,0)$ ，那么该点的梯度角度将对这四个子区域都有贡献。同时，假设该点的梯度角度（旋转以后）将落在第 $k=\lfloor 4\theta/\pi\rfloor$ 个bin，则该点的梯度角度对histogram中的第 $k$ 个和第 $k+1$ 个bin都将有贡献。也就是说一个像素点，最多会对描述子的8个值有影响。

subregionContribution.png

每个邻域内的像素点对histogram值的具体贡献值由梯度幅值加权而来，权重来自于几个部分。第一个部分是描述该位置与特征点相对位置关系的，离得越近权重越高，对应是一个高斯分布，其标准差为2，相对中心的位置是距离除以 $3\sigma$
$w_1 = \exp\left({-\frac{\left(\frac{x'}{3\sigma}\right)^2 + \left(\frac{y'}{3\sigma}\right)^2}{8}}\right)$
第二个部分则是描述该位置与该子区域中心位置的相对位置关系，离得越近权重越高。

weight.png

我们将上图中的蓝色阴影部分放大，阴影的四个角点对应了受影响的4个子区域的中心。我们定义该像素点到需要计算贡献的子区域的中心归一化距离分别为 $d_x,d_y$ ，则第二部分的权重为
$w_2 = d_x*d_y$
实际计算时，我们只需要计算到其中一个中心点的归一化距离，其他中心点的归一化距离可以由上图直接给出。

第三个部分是描述该点的梯度角度与bin的上下限的相对大小：
$\begin{eqnarray*} w_{3,k} &=& \frac{4\theta}{\pi} - \lfloor\frac{4\theta}{\pi}\rfloor\\ w_{3,k+1} &=& 1-\frac{4\theta}{\pi} + \lfloor\frac{4\theta}{\pi}\rfloor \end{eqnarray*}$

最后，某个像素点对histogram的贡献就可以写成
$h = m_G(x,y)*w_1*w_2*w_3$
最终得到的特征向量可以由下图表述

descriptor.png

描述子的归一化

为了应对不同亮度下的图片的差异，Sift特征会针对描述子向量进行一次 $L_2$ 归一化，这就是Sift不受亮度影响的原因。归一化之后，会根据阈值将描述向量里的一些低于阈值的数置成0。由于在上述加权中，高斯部分 $w_1$ 并没有归一化，所以有必要将归一化的 $L_2$ 范数与参与计算的像素点数相乘以后再与该阈值比较。在VLFeat中，该阈值设成了0，即这部分是disable的。

另外，还需要对描述子向量的值做一个极大值的抑制，VLFeat中设置的阈值为0.2，即对于大于0.2的描述子向量分量，直接将其设成0.2。由于极大值的抑制，描述子向量的范数将有可能不再是1，于是，我们需要对描述子再做一次 $L_2$ 归一化。

比较神奇的是，在Colmap的实现中，对上述已经 $L_2$ 归一化后的描述子向量还做了一次 $L_1-$ Root归一化
$l_i = \frac{h_i}{\sqrt{\sum_i h_i}}$
在Colmap中，最后会将float的描述子向量转化成unsigned byte，具体而言就是乘以512再取整。

DSP-Sift特征

DSP是domain size pooling的缩写。它与标准的Sift特征的区别在于描述子，即它只改变了特征点描述的方式，并不会改变特征点的生成。在描述子向量的生成过程中，我们设定的邻域范围是每个子区域均为 $3\sigma$ ，而对于DSP而言，它会在空间再进行一次采样。即它会将邻域范围乘上一个scale，在VLFeat的实现中，这个scale从1/3（对应子区域大小为 $\sigma$ ）到6（对应子区域大小为 $18\sigma$ ），中间10等分。对于每个scale，可以按照上述方法计算一个描述子向量，最后将上述10个描述子向量直接平均，再经过上述归一化的过程，即得到了DSP-Sift描述子。在我们的实验中，DSP-Sift确实能提高一些稀疏重建的匹配性能。

最后编辑于：2021.12.29 19:50:54

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 212,332评论 6赞 493
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 90,508评论 3赞 385
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 157,812评论 0赞 348
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 56,607评论 1赞 284
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 65,728评论 6赞 386
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 49,919评论 1赞 290
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,071评论 3赞 410
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 37,802评论 0赞 268
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,256评论 1赞 303
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 36,576评论 2赞 327
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 38,712评论 1赞 341
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 34,389评论 4赞 332
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,032评论 3赞 316
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,798评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,026评论 1赞 266
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 46,473评论 2赞 360
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 43,606评论 2赞 350

SIFT特征详解

特征点的确定

高斯模糊

金字塔的构建

Octave数目的确定

Sublevel 数目的确定

每一个图像对应的模糊尺度

高斯差分图

极值的寻找

极值的调整

极值的亚像素定位

极值的筛选

特征点的尺度

确定特征点方向

梯度图的选取

邻域大小的确定

角度方向直方图的生成

主方向的确立

特征点的输出

特征点描述子

邻域大小的确定

角度方向直方图的计算

描述子的归一化

DSP-Sift特征

推荐阅读更多精彩内容

每一个图像对应的模糊尺度 $\sigma$