RobHess的SIFT代码解析步骤一

平台：win10 x64 +VS 2015专业版 +opencv-2.4.11 + gtk_-bundle_2.24.10_win32

主要参考：1.代码：RobHess的SIFT源码：SIFT+KD树+BBF算法+RANSAC算法

2.书：王永明 王贵锦 《图像局部不变性特征与描述》

 

研究代码之前的几个问题issue：

1.问题描述：RobHess和David.Lowe关于SIFT的关系是什么，David.Lowe初创了SIFT算法，但是网络上为什么流传的更多的是RobHess的SIFT源码？为什么研究RobHess的SIFT源码？

答：参考期刊：王冬，夏乙等.基于SIFT特征的图像匹配算法实现.计算机与数字工程.2013(281)

 

2.问题描述：如何实现RobHess的SIFT源码？去哪里下载？怎么配置？需要下载哪些软件？

答案：RobHess的SIFT源码下载地址：http://robwhess.github.io/opensift/

配置及软件下载具体参看我的另一篇博客：结合OPENSIFT源码详解SIFT算法

 

3.问题描述：RobHess的SIFT源码解析有哪些比较好的博客或资源可以参考？

答案：（1）图像特征之SIFT：https://blog.csdn.net/u014568921/article/details/52514181

推荐原因：总述类型的（里边包括好多链接：综述；imgfeatures.h和imgfeatures.c文件；kdtree.h和kdtree.c文件；sift.h和sift.c文件；xform.h和xform.c文件）

（2）RobHess的SIFT源码分析：综述：https://blog.csdn.net/masibuaa/article/details/9191309

推荐原因：给了RobHess的SIFT源码，介绍了RobHess的SIFT源码中的几个文件都是做什么用的，方便刚接触的同学学习；SIFT算法的原理；k-d树算法的讲解；RANSAC算法。

（3）九之再续：教你一步一步用c语言实现sift算法、上：https://blog.csdn.net/v_JULY_v/article/details/6245939

推荐原因：在他写的关于sift算法的前倆篇文章里头，已经对sift算法有了初步的介绍：九、图像特征提取与匹配之SIFT算法，而后在：九（续）、sift算法的编译与实现里，也简单记录下了如何利用opencv，gsl等库编译运行sift程序。

缺点：博客作者代码已经缺失，但推荐加群：169056165，里边有相关的群资料

（4）【特征匹配】SIFT原理与C源代码剖析：https://www.cnblogs.com/liguangsunls/p/6785878.html

推荐原因：构建特征描写叙述子这块的公式（搜索半径，旋转公式，权值累计）比较详细。

（5）SIFT算法实现及代码详解：https://wenku.baidu.com/view/11a564da5022aaea998f0f1e.html

 推荐原因：给出了存储在文件（feature1.txt）中的SIFT特征分析。

 

 

4.问题描述：RobHess的SIFT源码如何复现？个人电脑环境如何配置？

答案：参看我的博客：RobHess的SIFT环境配置：https://www.cnblogs.com/Alliswell-WP/p/RobHessSIFTOfEnvironment.html

 

5.问题描述：SIFT开源代码有哪些？

答案：
1）http://www.vlfeat.org vlfeat开源图像处理库，其中http://www.vlfeat.org/api/sift.html关于其代码中一些细节和SIFT原理的解释。
2）http://robwhess.github.io/opensift/

RobHess sift
3）http://www.cs.ubc.ca/~lowe/research.html

  David Lowe Research Projects中的SIFT
4）http://docs.opencv.org/2.4/modules/nonfree/doc/feature_detection.html

OpenCV2.4以后的版本已实现了SIFT，其源码和RobHess的很相似。

 

研究代码中的问题issue汇总：

一、RobHess的SIFT源码中的几个文件说明？

(1) imgfeatures.h和imgfeatures.c文件
imgfeatures.h中有SIFT特征点结构struct feature的定义，这个结构很重要，后面都要用到，除此之外还有一些特征点的导入导出以及特征点绘制函数的声明。
对应的imgfeatures.c文件中是特征点的导入导出以及特征点绘制函数的实现。

(2) utils.h和utils.c文件
这两个文件中是一些图像基本操作的函数，包括：
1、获取某位置的像素点
2、设置某位置的像素点（8位，32位和64位），
3、计算两点之间的距离的平方
4、在图片某一点画一个“X”
5、将两张图片合成为一个(在特征匹配中用到)，高是二者之和，宽是二者的较大者。

(3) sift.h和sift.c文件
这两个是最重要的，里面的内容说白了很简单，就是两个特征点检测函数sift_features()和 _sift_features()，
sift_features()是用默认参数进行特征点检测， _sift_features()允许用户输入各种检测参数，其实sift_features()中也是再次调用_sift_features()函数。
所以，你只需提供原图像和存储特征点的数组以及其他一些检测参数，然后调用sift_features()或 _sift_features()就可完成SIFT特征点检测。

但是这两个文件分析起来也是最复杂的。
sift.h中有默认的各种特征检测中用到的参数的宏定义，sift.c中是检测函数的实现，是最核心的一个文件，里面的一些未暴露接口的本地函数就有31个之多。
仔细分析完后发现大神就是大神，程序写的条理非常清楚，尤其是函数_sift_features()的函数体，从头到尾的几部分正好对应SIFT算法的几个步骤，把一些细节全部放在子函数中，非常清楚明了。

(4) minpq.h和minpq.c文件
这两个文件中实现了最小优先级队列(Minimizing Priority Queue)，也就是小顶堆，在k-d树的建立和搜索过程中要用到。

(5) kdtree.h和kdtree.c文件
这两个文件中实现了k-d树的建立以及用BBF(Best Bin First)算法搜索匹配点的函数。
如果你需要对两个图片中的特征点进行匹配，就要用到这两个文件。

(6) xform.h和xform.c文件
这两个文件中实现了RANSAC算法(RANdom SAmple Consensus 随机抽样一致)。
RANSAC算法可用来筛选两个图像间的SIFT特征匹配并计算变换矩阵。

(7) 其他文件：dspfeat.c，match.c，siftfeat.c
这几个文件中是一些使用SIFT算法的小例子：
dspfeat.c文件可以从预先保存的特征点文件中读取特征点并显示在图片上。
match.c文件可以检测两个图像中的特征点并进行匹配。
siftfeat.c文件利用SIFT算法检测特征点，有一些控制台操作提示。

 

二、RobHess的SIFT源码中的宏定义说明？

(1) imgfeatures.h和imgfeatures.c文件
/*画出特征点的颜色 */
#define FEATURE_OXFD_COLOR CV_RGB(255,255,0)
#define FEATURE_LOWE_COLOR CV_RGB(255,0,255)

#define FEATURE_MAX_D 128 //最大特征描述子的长度，定为128
 
//圆周率
#define M_PI        3.14159265358979323846

(2) utils.h和utils.c文件
#define ABS(x) ( ( (x) < 0 )? -(x) : (x) ) //绝对值函数

(3) sift.h和sift.c文件
//高斯金字塔每组内的层数，S=3，书P83
#define SIFT_INTVLS 3
 
//第0层的初始尺度，即第0层高斯模糊所使用的参数，σo=1.6，书P83
#define SIFT_SIGMA 1.6
 
//对比度阈值 |D(x)|，针对归一化后的图像，用来去除不稳定特征，|D(x)|<0.03，书P87
#define SIFT_CONTR_THR 0.04 //此处与书上的不太一样，RobHess此处的阈值使用的是0.04/S
 
//主曲率比值的阈值，用来去除边缘特征，r=10，书P88
#define SIFT_CURV_THR 10
 
//是否将图像放大为之前的两倍，书P82，图像在计算高斯尺度空间前先将尺度扩大一倍（第0组o=-1那层索引）
#define SIFT_IMG_DBL 1
 
 
//计算特征描述子过程中，计算方向直方图时，将特征点附近划分为d*d个区域，每个区域生成一个直方图，SIFT_DESCR_WIDTH即d的默认值，描述子直方图矩阵默认尺寸Bp=4，书P133
#define SIFT_DESCR_WIDTH 4
 
//计算特征描述子过程中，每个方向直方图的bin个数，8个方向的梯度直方图，书P133
#define SIFT_DESCR_HIST_BINS 8

//输入图像的尺度为0.5，高斯模糊/掩模，书P81
#define SIFT_INIT_SIGMA 0.5
 
//边界的像素宽度，检测过程中将忽略边界线中的极值点，即只检测边界线以内是否存在极值点
#define SIFT_IMG_BORDER 5 //书中未涉及
 
//通过插值进行极值点精确定位时，最大插值次数，即关键点修正次数
#define SIFT_MAX_INTERP_STEPS 5 //书中未涉及最多插值多少次
 
//特征点方向赋值过程中，梯度方向直方图中柱子(bin)的个数，0度~360度分36个柱，每10度一个柱，书P130
#define SIFT_ORI_HIST_BINS 36
 
//特征点方向赋值过程中，搜索邻域的半径为：3 * 1.5 * σ，分配方向时高斯权重参数因子，高斯加权（1.5σ），书P132
#define SIFT_ORI_SIG_FCTR 1.5
 
//特征点方向赋值过程中，搜索邻域的半径为：3 * 1.5 * σ，分配方向时计算区域大小，根据3sigma外可忽略，圆半径3σ（3*1.5σ）图6-2，书P131
#define SIFT_ORI_RADIUS 3.0 * SIFT_ORI_SIG_FCTR
 
//特征点方向赋值过程中，梯度方向直方图的平滑次数，计算出梯度直方图后还要进行高斯平滑
#define SIFT_ORI_SMOOTH_PASSES 2 //书中未涉及具体多少次平滑处理
 
//特征点方向赋值过程中，梯度幅值达到最大值的80%则分裂为两个特征点，相当于主峰值80%能量的峰值，书P132
#define SIFT_ORI_PEAK_RATIO 0.8

//计算特征描述子过程中，特征点周围的d*d个区域中，每个区域的宽度为m*σ个像素，SIFT_DESCR_SCL_FCTR即m的默认值，σ为特征点的尺度，m=3，书P133
#define SIFT_DESCR_SCL_FCTR 3.0
 
//计算特征描述子过程中，特征描述子向量中元素的阈值(最大值，并且是针对归一化后的特征描述子)，超过此阈值的元素被强行赋值为此阈值，归一化处理中，对特征矢量中值大于0.2进行截断（大于0.2只取0.2），书P135
#define SIFT_DESCR_MAG_THR 0.2
 
//计算特征描述子过程中，将浮点型的特征描述子变为整型时乘以的系数
#define SIFT_INT_DESCR_FCTR 512.0 //书中未涉及
 
//定义了一个带参数的函数宏，用来提取参数f中的feature_data成员并转换为detection_data格式的指针，返回特征检测数据
#define feat_detection_data(f) ( (struct detection_data*)(f->feature_data) )


(4) minpq.h和minpq.c文件
//要为其分配空间的优先级队列元素的初始值
#define MINPQ_INIT_NALLOCD 512

(5) kdtree.h和kdtree.c文件
无

(6) xform.h和xform.c文件
/*RANSAC算法的容错度
对于匹配点对<pt,mpt>，以及变换矩阵H，
如果pt经H变换后的点和mpt之间的距离的平方小于RANSAC_ERR_TOL，则可将其加入当前一致集
*/
#define RANSAC_ERR_TOL 3
 
//内点数目占样本总数目的百分比的最小值
#define RANSAC_INLIER_FRAC_EST 0.25
 
//一个匹配点对支持错误模型的概率（不知道是干什么用的）
#define RANSAC_PROB_BAD_SUPP 0.10
 
//定义了一个带参数的函数宏，用来提取参数feat中的feature_data成员并转换为ransac_data格式的指针
#define feat_ransac_data( feat ) ( (struct ransac_data*) (feat)->feature_data )
 

(7) match.c文件

// 在k-d树上进行BBF搜索的最大次数，kdtree_bbf_knn函数用到

#define KDTREE_BBF_MAX_NN_CHKS 200  //书P158超时检查机制，具体多少，书中未涉及

 

/*  最近点与次最近点距离之比*/

//目标点与最近邻和次近邻的距离的比值的阈值，若大于此阈值，则剔除此匹配点对

//通常此值取0.6，值越小找到的匹配点对越精确，但匹配数目越少

#define NN_SQ_DIST_RATIO_THR 0.49  //与代码不同，判断阈值THR=0.49，书P163

 

三、RobHess的SIFT源码中使用的结构体及存储说明？按用到的先后顺序排列

（1）imgfeatures.h——结构体feature；feature_type；feature_match_type；
　　 imgfeatures.c——无

/*特征点结构体
此结构体可存储两种类型的特征点：
FEATURE_OXFD表示是牛津大学VGG提供的源码中的特征点格式，
FEATURE_LOWE表示是David.Lowe提供的源码中的特征点格式。
如果是OXFD类型的特征点，结构体中的a,b,c成员描述了特征点周围的仿射区域(椭圆的参数)，即邻域。
如果是LOWE类型的特征点，结构体中的scl和ori成员描述了特征点的大小和方向。
fwd_match，bck_match，mdl_match一般同时只有一个起作用，用来指明此特征点对应的匹配点
*/

struct feature

{

  double x;                 /*特征点x的坐标 */

  double y;                  /*特征点y的坐标*/

  double a;                 /* OXFD特征点中椭圆的参数 */

  double b;                 /*OXFD特征点中椭圆的参数*/

  double c;                /*OXFD特征点中椭圆的参数 */

  double scl;               /*LOWE特征点的尺度*/

  double ori;               /*LOWE特征点的方向*/

  int d;               /*特征描述子的长度，即维数，一般是128*/

  double descr[FEATURE_MAX_D];   /*128维的特征描述子，即一个double数组*/

  int type;              /*特征点类型 OXFD or LOWE */

  int category;            /*通用功能类别 */

  struct feature* fwd_match;     /*指明此特征点对应的匹配点*/

  struct feature* bck_match;    /*指明此特征点对应的匹配点*/

  struct feature* mdl_match;    /*指明此特征点对应的匹配点fwd_match，bck_match，mdl_match一般同时只有一个起作用，用来指明此特征点对应的匹配点*/

  CvPoint2D64f img_pt;       /*特征点的坐标,等于(x,y)*/

  CvPoint2D64f mdl_pt;       /*当匹配类型是mdl_match时用到*/

  void* feature_data;       /*用户定义的数据:
                  //在SIFT极值点检测中，是detection_data结构的指针
                 //在k-d树搜索中，是bbf_data结构的指针
                 //在RANSAC算法中，是ransac_data结构的指针*/

};

 

/*特征点的类型：
FEATURE_OXFD表示是牛津大学VGG提供的源码中的特征点格式，
FEATURE_LOWE表示是David.Lowe提供的源码中的特征点格式
*/
/** FEATURE_OXFD <BR> FEATURE_LOWE */
enum feature_type
{
    FEATURE_OXFD,
    FEATURE_LOWE,
};

 

/*特征点匹配类型：
FEATURE_FWD_MATCH：表明feature结构中的fwd_match域是对应的匹配点
FEATURE_BCK_MATCH：表明feature结构中的bck_match域是对应的匹配点
FEATURE_MDL_MATCH：表明feature结构中的mdl_match域是对应的匹配点
*/
enum feature_match_type
{
    FEATURE_FWD_MATCH,
    FEATURE_BCK_MATCH,
    FEATURE_MDL_MATCH,
};

 

 

（2）utils.h——无
　　 utils.c——无

 

（3）sift.h——结构体detection_data
　　 sift.c——无

struct detection_data //保存与检测有关的特征数据

{

  int r; //特征点所在的行

  int c; //特征点所在的列

  int octv; //特征点所在的组

  int intvl; //特征点所在的组中的层

  double subintvl; //特征点层方向（sigma方向，即intval方向）的亚像素偏移量

  double scl_octv; //特征点所在的组的尺度

};

 

（4）minpq.h——结构体min_pq和pq_node
　　 minpq.c——无

struct pq_node //最小化优先级队列中的元素

{

  void* data;

  int key;

};

struct min_pq //最小化优先级队列

{

  struct pq_node* pq_array;    /*数组包含优先级队列 */

  int nallocd;                 /* 分配的元素数 */

  int n;                       /*pq中的元素数量*/

};

 

（5）kdtree.h——结构体kd_node
　　 kdtree.c文件——结构体bbf_data

struct kd_node //K-D树中的结点结构

{

  int ki;                      /*分割位置(枢轴)的维数索引(哪一维是分割位置)，取值为1-128*/

  double kv;                   /*枢轴的值(所有特征向量在枢轴索引维数上的分量的中值)*/

  int leaf;                    /*是否叶子结点的标志 1 if node is a leaf, 0 otherwise */

  struct feature* features;    /*此结点对应的特征点集合(数组)*/

  int n;                       /*特征点的个数*/

  struct kd_node* kd_left;     /*左子树*/

  struct kd_node* kd_right;    /*右子树*/

};

 

struct bbf_data

{

  double d;

  void* old_data;

};

 

（6）xform.h——结构体ransac_data
　　 xform.c——无

//RANSAC算法中用到的结构

//在RANSAC算法过程中，此类型数据会被赋值给feature结构的feature_data成员

struct ransac_data

{

void* orig_feat_data; //保存此特征点的feature_data域的以前的值

int sampled; //标识位，值为1标识此特征点是否被选为样本

};

 

四、RobHess的SIFT源码中使用的IplImage结构体说明？

 参考：IplImage结构体介绍：https://blog.csdn.net/gflytu/article/details/46563159

IplImage结构体：https://wenku.baidu.com/view/cc937992a417866fb94a8ea0.html

 

SIFT四步骤和特征匹配及筛选：

步骤一：建立尺度空间，即建立高斯差分(DoG)金字塔dog_pyr

步骤二：在尺度空间中检测极值点，并进行精确定位和筛选创建默认大小的内存存储器

步骤三：特征点方向赋值，完成此步骤后，每个特征点有三个信息：位置、尺度、方向

步骤四：计算特征描述子

SIFT后特征匹配：KD树+BBF算法

SIFT后特征匹配后错误点筛选：RANSAC算法

 

 

步骤一：建立尺度空间，即建立高斯差分(DoG)金字塔dog_pyr

目录：
1、建立高斯尺度空间金字塔（GSS - Gauss Scale Space）
2、建立高斯差分金字塔（DOG - Difference of Gauss）

 

1.建立高斯尺度空间金字塔（GSS - Gauss Scale Space）

（1）问题描述：为什么使用高斯尺度空间去构造尺度不变性？怎么考虑的？

答：Koenderink（1984）和Lindeberg（1994）已经证明，在各种合理的假设下，唯一能产生尺度空间的核为高斯核函数，所以我们将图像的尺度空间表示成一个函数L(x,y,σ)，它是由一个变尺度的高斯函数G(x,y,σ)与图像I(x,y)卷积产生的。即L(x, y, σ) = G(x, y, σ) ∗ I(x, y)

其中*代表卷积操作，G(x, y, σ) 为二维高斯核函数，表示为：

 

（2）问题描述：目前的尺度分析的方法有哪些？只有图像金字塔吗？

答：参考硕士论文第10页：基于学习的图像多尺度表达的研究_阳平（北京工业大学）

可见，目前的多尺度分析方法主要包括图像金字塔结构方法、基于小波变换的多尺度分析方法以及多尺度几何分析。

 

（3）问题描述：如何构建高斯金字塔，代码如何解读？

答：代码描述：

/*将原图转换为32位灰度图并归一化，然后进行一次高斯平滑，

    并根据参数img_dbl决定是否将图像尺寸放大为原图的2倍*/

    init_img = create_init_img( img, img_dbl, sigma );

    //计算高斯金字塔的组数octvs，

    octvs = log( MIN( init_img->width, init_img->height ) ) / log(2) - 2;

    //为了保证连续性，在每一层的顶层继续用高斯模糊生成3幅图像，所以高斯金字塔每组

    //有intvls+3层，DOG金字塔每组有intvls+2层

    //建立高斯金字塔gauss_pyr，是一个octvs*(intvls+3)的图像数组

    gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );

分析：3.1）为什么需要图像需要放大为原图的2倍？怎么放大？

答案：LOWE建议把初始图像放大二倍，能够得到很多其它细节，而且觉得图像在SIFT_INIT_SIGMA=0.5时最清晰，初始高斯尺度为sigm=1.6。

放大步骤：

//调用函数，将输入图像转换为32位灰度图，并归一化

gray = convert_to_gray32( img );

#我理解：将3通道8位彩色图或8位灰度图转换为32位灰度图进行处理，因为32位颜色分辨率大大高于8位单通道的

//将图像长宽扩展一倍时，便有了底-1层，该层尺度为：

sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4 );

##sig_diff = sqrt( 1.6 * 1.6 - 0.5 *0.5 * 4 )=1.25

//创建放大图像，调用的opencv的函数cvCreateImage函数表示创建图像首地址并分配存储空间，cvSize函数表示矩形尺寸的结构，IPL_DEPTH_32F为32位单精度浮点数，1为单通道

dbl = cvCreateImage( cvSize( img->width*2, img->height*2 ),IPL_DEPTH_32F, 1 );

//放大原图的尺寸，gray为输入图像，dbl为输出图像，CV_INTER_CUBIC为插值方法——立方插值

cvResize( gray, dbl, CV_INTER_CUBIC );

//高斯平滑，调用的opencv的cvSmooth函数，高斯核在x,y方向上的标准差都是sig_diff，dbl为输入图像，dbl为输出图像，CV_GAUSSIAN为平滑方法——高斯平滑

cvSmooth( dbl, dbl, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );

 

3.2）convert_to_gray32( IplImage* img )如何将输入图像转换为32位灰度图,并进行归一化？

/*将输入图像转换为32位灰度图,并进行归一化

参数：

img：输入图像，3通道8位彩色图(BGR)或8位灰度图

返回值：32位灰度图*/

static IplImage* convert_to_gray32( IplImage* img )

{

    IplImage* gray8, * gray32;

 

gray32 = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

//首先将原图转换为8位单通道图像

if( img->nChannels == 1 )//若原图本身就是单通道，直接克隆原图，调用的opencv的cvClone函数

   gray8 = cvClone( img );

else//若原图是3通道图像（彩色图像）

{

    gray8 = cvCreateImage( cvGetSize(img), IPL_DEPTH_8U, 1 );//创建8位单通道图像

     cvCvtColor( img, gray8, CV_BGR2GRAY );//调用的opencv的cvCvtColor函数，将原图转换为8位单通道图像

}

//调用的opencv的cvConvertScale函数，然后将8位单通道图像gray8转换为32位单通道图像，并进行归一化处理(除以255)

cvConvertScale( gray8, gray32, 1.0 / 255.0, 0 );

 

//释放临时图像

cvReleaseImage( &gray8 );

 //返回32位单通道图像

return gray32;

}

 

3.3）怎么计算高斯金字塔的组数octvs？

答：

当输入850*680分辨率的图片，经过扩大一倍后为1700*1360，min(1700,1360)=1360，log(1360)/log(2)约等于10.4,10.4-2=8.4，因为octvs为int类型，直接去尾法等于8。

参考：https://blog.csdn.net/lingyunxianhe/article/details/79063547中1.1、高斯金字塔

 

3.4）问题描述：cvSmooth前两个参数被设置为0,窗口尺寸是多少？

答：函数原型说明：

void cvSmooth( const CvArr* src, CvArr* dst,
               int smoothtype=CV_GAUSSIAN,
               int param1=3, int param2=0, double param3=0 );
src
    输入图像.
dst
    输出图像.
smoothtype
    平滑方法:
        CV_BLUR_NO_SCALE (简单不带尺度变换的模糊) - 对每个象素领域 param1×param2 求和。如果邻域大小是变化的，可以事先利用函数 cvIntegral 计算积分图像。
        CV_BLUR (simple blur) - 对每个象素邻域 param1×param2 求和并做尺度变换 1/(param1•param2).
        CV_GAUSSIAN (gaussian blur) - 对图像进行核大小为 param1×param2 的高斯卷积
        CV_MEDIAN (median blur) - 发现邻域 param1×param1 的中值 (i.e. 邻域是方的).
        CV_BILATERAL (双滤波) - 应用双向 3x3 滤波，彩色 sigma=param1，空间 sigma=param2. 关于双向滤波，可参考http://www.dai.ed.ac.uk/CVonline/LOCAL_COPIES/MANDUCHI1/Bilateral_Filtering.html
param1
    平滑操作的第一个参数.代表滤波器窗口的宽度；
param2
    平滑操作的第二个参数. param2 为零对应简单的尺度变换和高斯模糊。代表滤波器窗口的高度；
param3
    对应高斯参数的 Gaussian sigma (标准差).代表水平方向的sigma值；

　 如果为零，这由下面的核尺寸计算：
    sigma = (n/2 - 1)*0.3 + 0.8, 其中 n=param1 对应水平核,
                                                    n=param2 对应垂直核.
    对小的卷积核 (3×3 to 7×7) 使用标准 sigma 速度会快。如果 param3 不为零，而 param1 和 param2 为零，则核大小有 sigma 计算 (以保证足够精确的操作).

param4
    对应高斯参数的 Gaussian sigma (标准差).代表垂直方向的sigma值；
注意：如果第三个，第四个参数已经指定，而前两个参数为0，那么窗口的尺寸由sigma确定（窗口的尺寸会根据sigma值自动确定），速度较慢但最有效

参看博客：1）https://blog.csdn.net/superjimmy/article/details/6179666

2）https://blog.csdn.net/wyl_steven/article/details/6333650

 

3.5.1）问题描述：gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );如何构建高斯金字塔

/*根据输入参数建立高斯金字塔

参数：

base：输入图像，作为高斯金字塔的基图像

octvs：高斯金字塔的组数

intvls：每组的层数

sigma：初始尺度

返回值：高斯金字塔，是一个octvs*(intvls+3)的图像数组*/

static IplImage*** build_gauss_pyr( IplImage* base, int octvs,

                                   int intvls, double sigma )

{

    IplImage*** gauss_pyr;

    //为了保证连续性，在每一层的顶层继续用高斯模糊生成3幅图像，所以高斯金字塔每组有

    //intvls+3层，DOG金字塔每组有intvls+2层

    double* sig = calloc( intvls + 3, sizeof(double));//每层的sigma参数数组

    double sig_total, sig_prev, k;

    int i, o;

 

    //为高斯金字塔gauss_pyr分配空间，共octvs个元素，每个元素是一组图像的首指针

　//calloc-C库函数在stdlib.h中，分配所需的内存空间，并返回一个指向它的指针，失败返回NULL

    gauss_pyr = calloc( octvs, sizeof( IplImage** ) );

    //为第i组图像gauss_pyr[i]分配空间，共intvls+3个元素，每个元素是一个图像指针

    for( i = 0; i < octvs; i++ )

        gauss_pyr[i] = calloc( intvls + 3, sizeof( IplImage* ) );

    //计算每次高斯模糊的sigma参数

    sig[0] = sigma;//初始尺度

    k = pow( 2.0, 1.0 / intvls );

    for( i = 1; i < intvls + 3; i++ )

    {

        sig_prev = pow( k, i - 1 ) * sigma;//i-1层的尺度

        sig_total = sig_prev * k;//i层的尺度

        sig[i] = sqrt( sig_total * sig_total - sig_prev * sig_prev );

    }

    //逐组逐层生成高斯金字塔

    for( o = 0; o < octvs; o++ )//遍历组

        for( i = 0; i < intvls + 3; i++ )//遍历层

        {

            if( o == 0  &&  i == 0 )//第0组，第0层，就是原图像

                gauss_pyr[o][i] = cvCloneImage(base);

            else if( i == 0 )//新的一组的首层图像是由上一组最后一层图像下采样得到

                gauss_pyr[o][i] = downsample( gauss_pyr[o-1][intvls] );

            else//对上一层图像进行高斯平滑得到当前层图像

            {   //创建图像

                gauss_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i-1]),IPL_DEPTH_32F, 1 );

                //对上一层图像gauss_pyr[o][i-1]进行参数为sig[i]的高斯平滑，得到当前层图像

                //gauss_pyr[o][i]

                cvSmooth( gauss_pyr[o][i-1], gauss_pyr[o][i], CV_GAUSSIAN, 0, 0, sig[i], sig[i] );

            }

        }

        free( sig );//释放sigma参数数组

        return gauss_pyr;//返回高斯金字塔

}

 

3.5.2）问题描述：每组尺度如何计算？每层尺度如何计算？我们在源代码上没有看到组间的2倍的关系？

1）所有空间尺度为： 

第1组2^0（σ，kσ，k^2 σ，……，k^(n-1) σ）=（σ，kσ，k^2 σ，……，k^(n-1) σ）　　　　  =（1.6，2^(1/3) *1.6，2^(2/3) *1.6，，2^(3/3) *1.6，2^(4/3) *1.6，2^(5/3) *1.6）

第2组2^1（σ，kσ，k^2 σ，……，k^(n-1) σ）=（2σ，2kσ，2k^2 σ，……，2k^(n-1) σ）　　 =（2*1.6，2*2^(1/3) *1.6，2*2^(2/3) *1.6，2*2^(3/3) *1.6，2*2^(4/3) *1.6，2*2^(5/3) *1.6）

第3组2^2（σ，kσ，k^2 σ，……，k^(n-1) σ）=（4σ，4kσ，4k^2 σ，……，4k^(n-1) σ）　　 =（4*1.6，4*2^(1/3) *1.6，4*2^(2/3) *1.6，4*2^(3/3) *1.6，4*2^(4/3) *1.6，4*2^(5/3) *1.6）

 

第4组2^3（σ，kσ，k^2 σ，……，k^(n-1) σ）=（8σ，8kσ，8k^2 σ，……，8k^(n-1) σ）　　 =（8*1.6，8*2^(1/3) *1.6，8*2^(2/3) *1.6，8*2^(3/3) *1.6，8*2^(4/3) *1.6，8*2^(5/3) *1.6）

第5组2^4（σ，kσ，k^2 σ，……，k^(n-1) σ）=（16σ，16kσ，16k^2 σ，……，16k^(n-1) σ）=（16*1.6，16*2^(1/3) *1.6，16*2^(2/3) *1.6，16*2^(3/3) *1.6，16*2^(4/3) *1.6，16*2^(5/3) *1.6）

 

第6组……

第7组……

第8组……

这个尺度因子都是在原图上进行的。而在算法实现过程中，采用高斯平滑是在上一层图像上再叠加高斯平滑。即我们在程序中看到的平滑因子为（4-5）

实际代码中的计算数值如我下边的演草纸上算的：

 

2）我们在源代码上同一时候也没有看到组间的2倍的关系，实际在对每组的平滑因子都是一样的，2倍的关系是因为在降採样的过程中产生的，第二层的第一张图是由第一层的平滑因子为2σ的图像上（即倒数第三张）降採样得到，此时图像平滑因子为σ，所以继续採用以上的平滑因子。

但在原图上看。形成了所有的空间尺度。

 

3）参考博客——尺度空间的连续性：https://blog.csdn.net/xw20084898/article/details/16832755

 

 3.6）问题描述： 3.5）构建高斯金字塔如何让下采样？

 答：代码描述及说明

/*对输入图像做下采样生成其四分之一大小的图像(每个维度上减半)，使用最近邻插值方法

参数：

img：输入图像

返回值：下采样后的图像*/

static IplImage* downsample( IplImage* img )

{

    //下采样图像

    IplImage* smaller = cvCreateImage( cvSize(img->width / 2, img->height / 2), img->depth, img->nChannels ); //创建比原图小四分之一大小的空矩阵(每个维度上减半)

    cvResize( img, smaller, CV_INTER_NN );//尺寸变换，调用的opencv的cvResize函数，img待处理图像, smaller处理后图像, CV_INTER_NN最近邻插值

    return smaller;

}

 用到的CV_INTER_NN最近邻插值如何实现？

 参看博客：

1）https://blog.csdn.net/zxj820/article/details/50594955

2）http://blog.chinaunix.net/uid-26020768-id-3187769.html

 

2、建立高斯差分金字塔（DOG - Difference of Gauss）

（1）问题描述：SIFT如何使用DoG？

答：SIFT算法建议，在某一尺度上对特征点的检测，可以通过对两个相邻高斯尺度空间的图像相减，得到一个DoG (Difference of Gaussians)的响应值图像D(x,y,σ)。然后，仿照LoG方法，通过对响应值图像D(x,y,σ)进行非最大值抑制(局部极大搜索，正最大和负最大)，在位置空间和尺度空间中定位特征点。其中

D(x, y, σ) = (G(x, y, kσ) − G(x, y, σ)) ∗ I(x, y) = L(x, y, kσ) − L(x, y, σ)

式中，k为相邻尺度空间倍数的常数。

 

（2）问题描述：为什么用DoG来检测特征点？原理是什么？

答：参考博客：https://www.cnblogs.com/JiePro/p/sift_1.html

Lindeberg证明用σ²标准化的高斯拉普拉斯(∇²G, LOG, Laplacian of Gauss)有着真正的尺度无关的特性，而Mikolajczyk发现，相比于其他一系列函数（比如梯度，Hessian，Harris角点函数等），用σ²标准化的高斯拉普拉斯(σ²∇²G)有着更稳定的图像特征，因此σ²∇²G函数是我们理想的寻找特征点的函数。而在此，我们利用DoG来近似σ²∇²G。

 

（3）问题描述： 为什么用前一个octave中的倒数第三幅图像生成下一octave中的第一幅图像？
答：实现了尺度的平滑过渡。详细参考博客——尺度空间的连续性：https://blog.csdn.net/xw20084898/article/details/16832755

（4）问题描述：dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );如何构建高斯差分(DoG)金字塔dog_pyr

答案：代码注释如下：

/*通过对高斯金字塔中每相邻两层图像相减来建立高斯差分金字塔

参数：

gauss_pyr：高斯金字塔

octvs：组数

intvls：每组的层数

返回值：高斯差分金字塔，是一个octvs*(intvls+2)的图像数组*/

static IplImage*** build_dog_pyr( IplImage*** gauss_pyr, int octvs, int intvls )

{

    IplImage*** dog_pyr;

    int i, o;

    //为高斯差分金字塔分配空间，共octvs个元素，每个元素是一组图像的首指针

    dog_pyr = calloc( octvs, sizeof( IplImage** ) );

    //为第i组图像dog_pyr[i]分配空间，共(intvls+2)个元素，每个元素是一个图像指针

    for( i = 0; i < octvs; i++ )

        dog_pyr[i] = calloc( intvls + 2, sizeof(IplImage*) );

    //逐组逐层计算差分图像

    for( o = 0; o < octvs; o++ )//遍历组

        for( i = 0; i < intvls + 2; i++ )//遍历层

        {   //创建DoG金字塔的第o组第i层的差分图像

            dog_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i]), IPL_DEPTH_32F, 1 );

            //将高斯金字塔的第o组第i+1层图像和第i层图像相减来得到DoG金字塔的第o组第i层图像

            cvSub( gauss_pyr[o][i+1], gauss_pyr[o][i], dog_pyr[o][i], NULL );//调用opencv的cvSub函数，gauss_pyr[o][i+1]被减矩阵, gauss_pyr[o][i]减矩阵, dog_pyr[o][i]结果矩阵

        }

        return dog_pyr;//返回高斯差分金字塔

}

 

（5）问题描述：高斯核性质及其在SIFT中的应用

答：参考：https://www.cnblogs.com/JiePro/p/sift_1.html

 对于二维高斯卷积，有如下性质：
1.距离高斯核中心3σ距离外的系数很小，相对于3σ内的系数值可以忽略不计，所以只用(6σ + 1)*(6σ + 1)的面积计算卷积即可。
2.线性可分，二维高斯核卷积（计算次数O(n2*M*N)）效果等于水平和竖直方向的两个一维高斯核（计算次数O(n*M*N) + O(n*M*N)））累积处理的效果。n*n为滤波器大小，M*N图像大小
3.对一幅图像进行多次连续高斯模糊的效果与一次更大的高斯模糊可以产生同样的效果，大的高斯模糊的半径是所用多个高斯模糊半径平方和的平方根。例如，使用半径分别为 6 和 8 的两次高斯模糊变换得到的效果等同于一次半径为 10 （勾股定理）的高斯模糊效果。根据这个关系，使用多个连续较小的高斯模糊处理不会比单个高斯较大处理时间要少。

其中，性质3有助我们更快速的生成GSS。虽然两个小半径处理时间并不会比单个高斯较大半径处理的时间少，但相对于原图像，我们已经有一个小半径模糊过的图像，即可用另一小半径在现成的图像上进行模糊得到较大半径的模糊效果。例如尺度为σ0的图像已经存在，我们要得到尺度为kσ0的图像，即可只使用尺度为[(kσ0)2-(σ0)2]1/2的核在尺度为σ0的图像上进行模糊，即可得到尺度为kσ0模糊的图像，同理可得到高斯尺度空间各层的图像，使得高斯尺度空间生产得更快。

 

 附：

（1）问题issue：malloc函数的一些问题？memset函数的一些问题？

参看：1)malloc函数百度百科：https://baike.baidu.com/item/malloc%E5%87%BD%E6%95%B0/8582146?fr=aladdin

 2)memset函数百度百科：https://baike.baidu.com/item/memset/4747579?fr=aladdin

 3)malloc、calloc、realloc的区别：https://www.cnblogs.com/lidabo/p/4611411.html

 

 

 

 

 

 

转载于:https://www.cnblogs.com/Alliswell-WP/p/SIFT_code1.html