SIFT算法计算原理

旧城等待， 2022-06-09 21:18 153阅读 0赞

特征提取在CV(computer vision)领域非常重要。SIFT是非常出名的特征提取算法，它来自论文IJCV'04的“Distinctive image features from scale-invariant keypoints”

本文中的SIFT实现在http://blogs.oregonstate.edu/hess/，http://blogs.oregonstate.edu/hess/code/sift/，Linux版本的。

## 0. 摘要 ##

这里说SIFT算法如何之好，在图像有尺度变化和旋转的时候，提取的特征是不变的，在一定程度的仿射扭曲，3D视角，噪声，照明变化的情况下，匹配稳定。

## 1. 介绍 ##

SIFT算法用级联式的过滤方法，降低提取特征的时间成本。也就是说，它用4个步骤提取特征，对图像上的每个像素，只有通过了着4个步骤，才算是特征点。这样做的好处就是，把大部分的像素在开始就排除掉，不在它们那里浪费时间。这4个步骤如下： **1）尺度（大小）空间*极值*；2）关键点*定位（坐标）*；3）方向*关联*；4）关键点*描述子*。**

## 2. 相关研究 ##

特征提取文献综述。

## （步骤1） ##

## 3. 尺度空间*极值*的检测 ##

开始进入正题了。尺度空间极值检测是4个步骤的第1个。为什么要检测尺度空间极值呢？如果要提取一个目标的不变特征，假设是我们要识别iphone，就要先给iphone拍张照片，然后用这个照片作为基准，提取特征。提取好特征值后，将来遇到iphone的大图，小图，旋转了一定角度的图，都能识别出来，如果只是大小和角度都不变的情况下识别，没什么实际价值。既然在大图，小图，旋转图的情况下都能识别，那么，必然在图像里必然存在一些“特性”，这些特征在大图，小图，旋转图里都“不变”的，或者近似不变。这个必然的，否则没有特征提取的必要了。尺度空间极值检测，就是找到图像上的若干候选点，这些点在图像的尺度发生变化的时候，也就是图像变大或者变小，它们和周围其他的点的关系保持不变。

引用的论文说，在一定假定下，尺度空间核只可能是高斯函数。公式（1）就给出了用微分高斯函数卷积图像的计算尺度空间极值的方法。注意：这块的算法本身是很简洁，只是写的不太容易看得懂，反正我看论文的话是没看明白如何实现它，第3节的位于3.1小节上的最后那段话写的相当费解。我从源代码出发，把第3节从开始到第3.1节之间的部分，把公式（1）的 **DOG计算方法** 写成如下 **步骤：**

\[1\] 假设待处理的图像是demo.jpg，是单通道的灰度图，尺寸是512\*256。

\[2\]  计算octave

什么是octave呢？octave（高斯金字塔每层）原来意思是音乐上的一个八度。代码里octave数量计算方法是：

octvs = log( MIN( init\_img->width, init\_img->height ) ) / log(2) - 2;

显然，如果图像是512\*256，按照上式计算octvs=6。

\[3\]  **计算图像高斯金字塔**

\[3.1\] 有 **两个参数intvls = 3, sigma = 1.6**，在代码里是常量。

\[3.2\] 计算高斯金字塔的函数，返回一个图像高斯金字塔的指针。**图像高斯金字塔指针**，指向了**6个octave指针**，每个**octave指针指向intvls+3=6个图像指针**。这也就是说，**图像指针是IplImage\***，那么，**每个octave指针就是指针的指针**，就是**IplImage \*\***，而 高斯金字塔的指针是指针的指针的指针，就是IplImage \*\*\* 。

\[3.3\] 计算sigma：在每个octave里，有intvls+3=6个图像，每个图像都是前一个图像通过高斯卷积获得的， **每次卷积的sigma不同**，需要分别计算 。sig\[0\] = sigma = 1.6。k = 2^(1/3)。其他的sig如下：

sig\[1\] = ( (k^0\*1.6\*k)^2 - ( k^0\*1.6 )^2   )^(0.5)

sig\[2\] = ( (k^1\*1.6\*k)^2 - ( k^1\*1.6 )^2   )^(0.5)

sig\[3\] = ( (k^2\*1.6\*k)^2 - ( k^2\*1.6 )^2   )^(0.5)

sig\[4\] = ( (k^3\*1.6\*k)^2 - ( k^3\*1.6 )^2   )^(0.5)

sig\[5\] = ( (k^4\*1.6\*k)^2 - ( k^4\*1.6 )^2   )^(0.5)

**对每个octave，sig都是相同的**。这个sig会用在每次octave的计算中。为什么有这么奇怪的计算方式呢？因为在源代码里，作者是根据同octave的第i个图像进行高斯卷积，计算第i+1个图像，才会有这么怪异的计算sigma，如果是用第1个计算其他的图像，就不会这么麻烦了。

\[3.4\] 计算过程

\[3.4.1\]第 0层 的octave第 0个 图像是 原始图像 demo.jpg。

\[3.4.2\]第 0层 的octave第 1个 图像是第0层octave的第 0个 图像的高斯卷积，卷积参数是sig\[1\]

\[3.4.3\]第 0层 的octave第2个图像是第0层octave的第1个图像的高斯卷积，卷积参数是sig\[2\]

\[3.4.4\]第0层的octave第3个图像是第0层octave的第2个图像的高斯卷积，卷积参数是sig\[3\]

\[3.4.5\]第0层的octave第4个图像是第0层octave的第3个图像的高斯卷积，卷积参数是sig\[4\]

\[3.4.6\]第0层的octave第 5个 图像是第0层octave的第 4个 图像的高斯卷积，卷积参数是sig\[5\]

\[3.4.7\]第 1层 的octave第 0个图 像是第 0层 octave的第 intvls=3个 图像的 下抽样，宽度和高度均原来一半 ，即256\*128。

\[3.4.8\]第 1层 的octave第 1个 图像是第 1层 octave的第**0个**图像的高斯卷积，卷积参数是sig\[1\]

... ...

以下以此类推，最终计算出 整**个高斯金字塔** 。

\[4\] 计算图像DOG金字塔

\[4.1\]DOG金字塔指针：DOG金字塔是一个指向指针的指针的指针，即IplImage \*\*\*。DOG金字塔指针，指向6个octave指针，每个octave指针，是指向指针的指针，即IplImage \*\*。每个octave指针指向5个图像指针。

\[4.2\] 计算过程：

\[4.2.1\]第 **0个**octave第 **0个**图像，是高斯金字塔第 **0个**octave第 **1个**图像 **减**其的第 **0个**octave第 **0个**图像（01-00）

\[4.2.2\]第0个octave第 **1个**图像，是高斯金字塔第 **0个**octave第 **2个**图像 **减**其的第 **0个**octave第 **1个**图像（02-01）

\[4.2.3\]第0个octave第 **2个**图像，是高斯金字塔第 **0个**octave第 **3个**图像 **减**其的第 **0个**octave第 **2个**图像（03-02）

\[4.2.4\]第0个octave第3个图像，是高斯金字塔第0个octave第4个图像减其的第0个octave第3个图像 （04-03）

\[4.2.5\]第0个octave第4个图像，是高斯金字塔第0个octave第5个图像减其的第0个octave第4个图像 （05-04）

\[4.2.6\]第1个octave第0个图像，是高斯金字塔第1个octave第1个图像减其的第1个octave第0个图像 （11-10）

\[4.2.7\]第1个octave第1个图像，是高斯金字塔第1个octave第2个图像减其的第1个octave第1个图像 （12-11）

\[4.2.8\]第1个octave第2个图像，是高斯金字塔第1个octave第3个图像减其的第1个octave第2个图像 （13-12）

\[4.2.9\]第1个octave第3个图像，是高斯金字塔第1个octave第4个图像减其的第1个octave第3个图像 （14-13）

... ...

![46090764_29.jpg][]

以下以此类推，最终计算出整个 **DOG金字塔** 。

### 3.1 **局部极值检测** ###

从DOG金字塔里检测尺度空间极值。 **DOG金字塔有6个octave（组数与高斯金字塔一样）**。对于每个octave，分别找极值。具体步骤如下：

\[5\] **计算尺度空间极值点**

\[5.1\] 第0个octave的尺度空间极值

\[5.1.1\]第 0层 空间极值：对于第 **1个** 图像上的 所有像素点的值 ，如果它比它的 8个邻域 和第 **0个** 图像的**9个邻域**像素和第 **2个** 图像的**9个邻域**像素都大，或者都小，那么它是一个极值点。

\[5.1.2\]第**1层**空间极值：对于第 **2个** 图像上的所有像素点的值，如果它比它的8个邻域和第 **1个** 图像的9个邻域像素和第 **3个** 图像的9个邻域像素都大，或者都小，那么它是一个极值点。

\[5.1.3\]第 2层 空间极值：对于第 **3个** 图像上的所有像素点的值，如果它比它的8个邻域和第 **2个** 图像的9个邻域像素和第 **4个** 图像的9个邻域像素都大，或者都小，那么它是一个极值点。

\[5.1.4\]第3层空间极值：对于第4个图像上的所有像素点的值，如果它比它的8个邻域和第3个图像的9个邻域像素和第5个图像的9个邻域像素都大，或者都小，那么它是一个极值点。

\[5.2\] 第 **1个** octave的尺度空间极值

... ...

以下以此类推，最终计算出所有的空间极值点。这些极值点是候选点，需要进一步 过滤 。

### 3.2 尺度抽样频率 ###

这一段分析抽样的各参数对性能的影响问题。

### 3.3 空间域的抽样频率 ###

这一段分析抽样的各参数对性能的影响问题。

## （步骤2） ##

## 4、 精确的关键点*定位* ##

这部分内容是 SIFT的第2步 。从第4节到第4.1节之间部分，是 精确定位关键的坐标 ，候选集里的所有点，都要进一步处理。每个关键点根据泰勒级数进行二阶展开，它的**像素之值是横坐标，纵坐标和sigma的函数**，要求它的极值点，就要将二阶泰勒展开，求一阶导数，令导数为0，求出det\_r, det\_c, det\_sigma，那么，就能 精确定位极值点 的精确的intvls，纵坐标和横坐标。求出精确坐标之后，再计算它的若干属性，如对比度，要超过一定的阈值才行。这就是公式（3）。对每个关键点，计算dD和hessen\_3D，然后矩阵相乘，得出det\_x，得出精确值。这一块对应的interp\_extremum函数，呆定的，熟悉线性代数和微积分的这块好理解，如果不好理解....翻书去看吧。注意，着一块的推导，类似ax^2+bx+c求一阶导数，将一阶偏导和二阶偏导视为常量即可，但如果将它们视为函数，结果是不同的，是(3)的二分之一。

\[6\]精确定位每个关键点的坐标

\[6.1\] 对每个关键点，用一个循环，计算关键点精确坐标，最多计算5次。

\[6.1.1\] 插值计算其关键点极值在intvsl, row, col上的增量xi, xr, xc。这里的计算是根据公式（3）进行的。

\[6.1.2\] 如果xi, xr, xc的绝对值都小于0.5，表明够精确了，跳出循环。

\[6.1.3\] 对intvl, r, c 进行更新。如果intvl, r, c超过正常值，就返回NULL。

\[6.2\] 如果循环超过给定值，返回NULL。

\[6.3\] 计算关键点的对比度，也就是第11页论文的公式。如果对比度小于某个常量，返回NULL。

\[6.4\] 生成关键点数据结构feat，根据octave的值，对关键点的坐标进行放大，计算其在原始图像上的像素位置。同时，也把关键点所在的octave, r, c, intvl, subintvl, 并返回feat。

4.1 删除边缘点。

如果关键点是边缘线条上的，那么要去除掉。 对所有的关键点都这么处理，一个公式，呆定的。

这个原理是这样的，如果一个点，它在线条边缘上，那么，它的2维hessian矩阵的特征值，必然是一个特征值比另外一个特征值大很多。如果它在 两个线条的交点，也就是角点 ，那么两个特征值在同一个数量级上。对hessian矩阵直接求特征值计算量比较大，因此，论文第12页用一个替代性的方法求解。

\[7\] 删除边缘上的关键点

对于所有的关键点，即\[6\]里头得到的每个feat，根据第12页的最后一个公式，判断是否惟边缘上的点，如果是，删除之。

\[8\] **计算特征属性尺度**

对\[7\]得到的每个特征，计算特征尺度，对应的函数是calc\_feature\_scales。在\[6.4\]里，要计算subintvl，这是一个小数，实际上，每个图像的octave和intvl都是整数，但整数的intvlb并不能精确表示极值点，subintvl就是极值点所在的intvl的增量，这一点在\[6.4\]步骤是很清楚的。在feat的数据结构里，有一个ddata数据结构，记录suvintvl，在这里用到它了。语句：feat->scl = sigma \* pow( 2.0, ddata->octv + intvl / intvls )中，signa是常量1.6，octv就是极值点所在的octave的序号，intvl是极值点所在的精确intvl，它很可能是带小数点的实数。ddata->scl\_octv = sigma \* pow( 2.0, intvl / intvls )就含义类似。这两个值在第3步，进行方向关联的时候，取关键点的邻域时候用得到。

##     ##

## （步骤1） ##

## **5. 方向*关联*** ##

每个关键点都跟关联一个常量的方向值，让特征描述子对图像旋转具有不变性。

\[9\] 计算特征方向

对应函数calc\_feature\_oris。

\[9.1\] 对于每个\[8\]中获取的特征，都要对它进行方向关联。也就是说，确定它所在的高斯金字塔图像的一定大小的邻域内的主方向。

\[9.2\]计算直方图函数是ori\_hist，输入参数是，该特征在高斯图像金字塔所在图像的指针，该特征在这个图像上的行号和列号，直方图的桶数(bins)，这里是常量36，该特征的邻域大小，注意，这里是根据\[8\]里计算出来的scl\_octv乘以一个常量计算的，也就是说，基本上，intvl线性变化，最后一个参数是sigma，也是根据intvl线性变化。这个函数的计算是呆定的，是按照论文上写的做的。需要注意的是：源代码使用的是atan2函数，计算出来的角度，跟纵横坐标的符号相关，它的值域在\[-pi, pi\]，这里做的处理是bin = cvRound( n \* ( ori + CV\_PI ) / PI2 )，ori+pi实在\[0, 2pi\]，也就是说，所有的角度都可能，每个bin里都会可能有值。计算出直方图后，对直方图进行平滑操作。

\[9.3\] 将“好的”特征压栈：对于直方图的每个bin，如果它比前一个bin和后一个bin都大且超过一个阈值，那么，就复制当前的feat，将这个bin对应的方向值跟feat做关联，然后压栈。注意，在步骤\[9.1\]是从栈的前面pop一个feat进行处理，而这里是将特征push，这种用法效率比开两个栈要高一些。此外，这种计算方法，一个特征点可能生成多个新特征，只要它的方向值足够大。

## （步骤1） ##

## **6.图像局部*描述子*** ##

这是SIFT算法的四个步骤的最后一个。它在源代码里对应函数compute\_descriptors。

###  ###

###             6.1计算局部描述子 ###

###             6.2计算描述直方图 ###

6.2.1本文将特征点的 近邻区分成4\*4个区块 ，在 每个区块里，有若干个像素点， 这些像素点都具有**不同的方向和幅值**，要将这些像素点的方向和幅值计算直方图，直方图的bin是8个，也就是说，将每个bin对应45度。

6.2.2对于每个特征点，首先，计算radius，也就是它的计算描述直方图的 **邻域的半径**。

6.2.3对于邻域内的每个点，要将这个点的当前坐标，转化到 **主方向**上去。这一块在代码里对应的是：

c\_rot = ( j \* cos\_t - i \* sin\_t ) / hist\_width;

r\_rot = ( j \* sin\_t + i \* cos\_t ) / hist\_width;

rbin = r\_rot + d / 2 - 0.5;

cbin = c\_rot + d / 2 - 0.5;

if( rbin > -1.0   &&   rbin < d   &&   cbin > -1.0   &&   cbin < d )

求c\_rot和r\_rot是比较费解的。从后面我们可以知道，r+i，表明 **i是像素纵坐标**上的变量，c+j，表明 **j是像素横坐标**上的变量。如果ori为0，那么，此时i,j不需要旋转。因为主方向不变。假设ori不是0，那么，就需要将现在的坐标系，旋转到ori的方向上去。假设有坐标为\[i, j\]， **r是极轴，sita是极角**，那么，在极坐标下，i,j 点在新坐标系下面，就是x\_new = r\*cos(sita+ori) = r\*cos(sita)cos(ori) - r\*sin(sita)sin(ori) = j\*cos(ori)-i\*sin(ori)，y\_new = r\*sin(sita+ori) = r\*sin(sita)cos(ori) + r\*cos(sita)sin(ori) = i\*cos(oir) + j\*sin(ori)。这个要根据图像坐标计算，也就是说，原点在左上角，x向右 为正，y向下为正。

为什么要除以hist\_width呢？在转换之后的新坐标系下，要看c\_rot和r\_rot跟hist\_width之间的关系。因为，计算区域是限制在4个hist\_width的，也就是说， **c\_rot和r\_row的值，应该在-2\*hist\_width和2\*hist\_width之间**。

为什么要计算rbin和cbin呢？我们用下面方式推理：if条件，其实就是说-1 < rbin < d <==> -1 < rbin<4 <==> -1 < r\_rot+2-0.5 < 4 <==> -2.5 < r\_rot < 2.5。而 r\_rot是新坐标值与hist\_width的比值。if判断， **其实就是判断r\_rot是不是在指定的范围内，也就是几个hist\_width。**

如果if判断通过了，就是计算该像素的方向和幅值。计算方向之后，要将方向减去ori，因为此时计算出来的是当前坐标系下的方向，不是新坐标系下的方向。然后，要根据新坐标系下的坐标，计算权重。

6.2.4计算每个局部区域的直方图，对应函数是interp\_hist\_entry。在这个函数，根据计算出的rbin和cbin，算出来这个点，应该属于4个\*4区块的那个区块，然后，再根据它的方向和幅值，将它累加到不同的bin里。这样，就得到了论文p15页图的右侧的小图，唯一对区别，就是这个图不是4\*4，是2\*2。

6.2.5将描述直方图转化描述子，对应函数hist\_to\_descr。这个就没什么好说啦，就是将4\*4\*8个值，逐次排列起来，然后做归一化，再转化成整数。

## 7. 目标识别的应用 ##

这个是见仁见智，不一而足，可以在原文之上做更多的发挥。

[46090764_29.jpg]: http://image79.360doc.com/DownloadImg/2014/10/1315/46090764_29.jpg