特征提取在CV(computer vision)领域非常重要。SIFT是非常出名的特征提取算法，它来自论文IJCV'04的“Distinctive image features from scale-invariant keypoints”，在scholar.google.com上查到的引用次数一万多次，很高了！我准备在这个帖子里，根据这篇论文和SIFT算法的一个开源实现，详细描述SIFT算法。
本文中的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. 相关研究
特征提取文献综述。
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个图像
[4.2.2]第0个octave第1个图像，是高斯金字塔第0个octave第2个图像减其的第0个octave第1个图像
[4.2.3]第0个octave第2个图像，是高斯金字塔第0个octave第3个图像减其的第0个octave第2个图像
[4.2.4]第0个octave第3个图像，是高斯金字塔第0个octave第4个图像减其的第0个octave第3个图像
[4.2.5]第0个octave第4个图像，是高斯金字塔第0个octave第5个图像减其的第0个octave第4个图像
[4.2.6]第1个octave第0个图像，是高斯金字塔第1个octave第1个图像减其的第1个octave第0个图像
[4.2.7]第1个octave第1个图像，是高斯金字塔第1个octave第2个图像减其的第1个octave第1个图像
[4.2.8]第1个octave第2个图像，是高斯金字塔第1个octave第3个图像减其的第1个octave第2个图像
[4.2.9]第1个octave第3个图像，是高斯金字塔第1个octave第4个图像减其的第1个octave第3个图像
... ...
以下以此类推，最终计算出整个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 空间域的抽样频率
这一段分析抽样的各参数对性能的影响问题。
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步，进行方向关联的时候，取关键点的邻域时候用得到。
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，这种用法效率比开两个栈要高一些。此外，这种计算方法，一个特征点可能生成多个新特征，只要它的方向值足够大。
6.图像局部描述子
这是SIFT算法的四个步骤的最后一个。它在源代码里对应函数compute_descriptors。
[10]计算局部描述子
[10.1计算描述直方图
[10.1.1]本文将特征点的近邻区分成4*4个区块，在每个区块里，有若干个像素点，这些像素点都具有不同的方向和幅值，要将这些像素点的方向和幅值计算直方图，直方图的bin是8个，也就是说，将每个bin对应45度。
[10.1.2]对于每个特征点，首先，计算radius，也就是它的计算描述直方图的邻域的半径。
[10.1.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，因为此时计算出来的是当前坐标系下的方向，不是新坐标系下的方向。然后，要根据新坐标系下的坐标，计算权重。
[10.1.4]计算每个局部区域的直方图，对应函数是interp_hist_entry。在这个函数，根据计算出的rbin和cbin，算出来这个点，应该属于4个*4区块的那个区块，然后，再根据它的方向和幅值，将它累加到不同的bin里。这样，就得到了论文p15页图的右侧的小图，唯一对区别，就是这个图不是4*4，是2*2。
[10.2]将描述直方图转化描述子，对应函数hist_to_descr。这个就没什么好说啦，就是将4*4*8个值，逐次排列起来，然后做归一化，再转化成整数。
7. 目标识别的应用
这个是见仁见智，不一而足，可以在原文之上做更多的发挥。