RobHess的SIFT代码解析步骤二
平台：win10 x64 +VS 2015专业版 +opencv-2.4.11 + gtk_-bundle_2.24.10_win32
主要参考：1.代码：RobHess的SIFT源码
2.书：王永明 王贵锦 《图像局部不变性特征与描述》
SIFT四步骤和特征匹配及筛选：
步骤一：建立尺度空间，即建立高斯差分(DoG)金字塔dog_pyr
步骤二：在尺度空间中检测极值点，并进行精确定位和筛选创建默认大小的内存存储器
步骤三：特征点方向赋值，完成此步骤后，每个特征点有三个信息：位置、尺度、方向
步骤二：在尺度空间中检测极值点，并进行精确定位和筛选创建默认大小的内存存储器
问题及解答：
（1）问题描述：如何找寻关键点？在几层之间对比？思路是什么？
答：极值的检测是在DoG空间进行的，检测是以前点为中心，3pixel*3pixel*3pixel的立方体为邻域，判断当前点是否为局部最大或最小。如下图所示，橘黄色为当前检测点，绿色点为其邻域。因为要比较当前点的上下层图像，所以极值检测从DoG每层的第2幅图像开始，终止于每层的倒数第2幅图像（第1幅没有下层，最后1幅没有上层，无法比较）。
参考书P84-P85及图4-6
（2）问题描述：为什么还要对找寻的关键点进行插值？为什么要精确定位？
答：以上极值点的搜索时在离散空间中进行的，检测到的极值点并不是真正意义上的极值点。如下图所示，连续空间中极值与离散空间的区别。通常通过插值的方式，利用离散的值来插值，求取接近真正的极值的点。
对于一维函数，利用泰勒级数，将其展开为二次函数：
f(x) ≈ f(0) + f'(0)x + f''(0)x2
对于二维函数，泰勒展开为：
矩阵表示为：
矩阵表示为：
矢量表示为：
当矢量为n维时，有：
书中P86页有错误，我已经更正，并附在：https://www.cnblogs.com/Alliswell-WP/p/SIFT.html
求取f(x)的极值，只需求取∂f/∂x = 0。对于极值，x,y,σ三个变量，即为三维空间。利用三维子像元插值，设其函数为D(x, y, σ)，令x = (x, y, σ)T，那么在第一节中找到的极值点进行泰勒展开为（式1）如下：
（式1）
其中D为极值点的值，∂DT/∂x为在极值点各自变量的倒数，∂2D/∂x2为其在展开点相应的矩阵。对上式求导，另∂D(x)/∂x = 0，结果如下式，对应的
（为了表示方便，ˆx代替
）向量即为真正极值点偏离插值点的量。求解得（式2）如下：
（式2）
最终极值点的位置即为插值点x+ˆx，且多次迭代可以提高精度（一般为5次迭代）。
参考书P85-P87及图4-7和图4-8
（3）问题描述：既然知道了检测到的极值点并不是真正意义上的极值点，那怎么怎么插值呢？
答：图像离散的像素点组成的，所以要差分代替偏导，计算一阶偏导，二阶偏导构建Hessian矩阵。
（4）问题描述：对关键点定位后，需要剔除一些不好的KeyPoint，那什么是不好的KeyPoint的呢？如何剔除？
答：
1.DoG响应较低的点，即极值较小的点。
2.响应较强的点也不是稳定的特征点。DoG对图像中的边缘有较强的响应值，所以落在图像边缘的点也不是稳定的特征点。一方面图像边缘上的点是很难定位的，具有定位的歧义性；另一方面这样的点很容易受到噪声的干扰变得不稳定。
对于第一种，只需计算矫正后的点的响应值D(ˆx)，响应值小于一定阈值，即认为该点效应较小，将其剔除。将（式2）带入（式1），求解得：
在Lowe文章中，将|D(ˆx)|<0.03（图像灰度归一化为[0,1]）的特征点剔除。
对于第二种，利用Hessian矩阵来剔除。一个平坦的DoG响应峰值在横跨边缘的地方有较大的主曲率，而在垂直边缘的地方有较小的主曲率。主曲率可以通过2×2的Hessian矩阵H求出：
D值可以通过求临近点差分得到。H的特征值与D的主曲率成正比，具体可参见Harris角点检测算法。为了避免求具体的值，我们可以通过H将特征值的比例表示出来。令
为最大特征值，
为最小特征值，那么：
Tr(H)表示矩阵H的迹，Det(H)表示H的行列式。令
表示最大特征值与最小特征值的比值，则有：
上式与两个特征值的比例有关。随着主曲率比值的增加，
也会增加。我们只需要去掉比率大于一定值的特征点。Lowe论文中去掉r=10的点。
参考书P87-P88
（5）问题描述：RobHess的SIFT源码如何实现步骤二，大概思路是这样的？
答：（5.1）代码及说明：
/*步骤二：在尺度空间中检测极值点，并进行精确定位和筛选创建默认大小的内存存储器*/
storage = cvCreateMemStorage( 0 ); //调用opencv的cvCreateMemStorage函数，函数功能：用来创建一个内存存储器，来统一管理各种动态对象的内存。函数返回一个新创建的内存存储器指针。
//在尺度空间中检测极值点，通过插值精确定位，去除低对比度的点，去除边缘点，
//返回检测到的特征点序列
features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr, curv_thr, storage );
//计算特征点序列features中每个特征点的尺度
calc_feature_scales( features, sigma, intvls );
if( img_dbl ) //若设置了将图像放大为原图的2倍
adjust_for_img_dbl( features );//将特征点序列中每个特征点的坐标减半
//(当设置了将图像放大为原图的2倍时，特征点检测完之后调用)
（5.2）scale_space_extrema代码及说明：
/*在尺度空间中检测极值点，通过插值（三维二次函数）精确定位，去除低对比度的点，去除边缘点，返回检测到的特征点序列
参数：
dog_pyr：高斯差分金字塔
octvs：高斯差分金字塔的组数
intvls：每组的层数
contr_thr：对比度阈值，针对归一化后的图像，用来去除不稳定特征
cur_thr：主曲率比值的阈值，用来去除边缘特征
storage：存储器
返回值：返回检测到的特征点的序列*/
static CvSeq* scale_space_extrema( IplImage*** dog_pyr, int octvs, int intvls,
double contr_thr, int curv_thr, CvMemStorage* storage ) //octvs=O(log(min(length,width))/log(2)-2)；intvls=3；curv_thr=0.04；curv_thr=10
{
CvSeq* features;//特征点序列
double prelim_contr_thr = 0.5 * contr_thr / intvls;//像素的对比度阈值，此处与LOWE论文不同，，|D(x)|<0.03，书P87
struct feature* feat;
struct detection_data* ddata;
int o, i, r, c;
//在存储器storage上创建存储极值点的序列，其中存储feature结构类型的数据
features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct feature), storage );
/*遍历高斯差分金字塔，检测极值点*/
//SIFT_IMG_BORDER指明边界宽度，只检测边界线以内的极值点
for( o = 0; o < octvs; o++ )//第o组
for( i = 1; i <= intvls; i++ )//遍i层
for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)//第r行
for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)//第c列
//进行初步的对比度检查，只有当归一化后的像素值大于对比度
//阈值prelim_contr_thr时才继续检测此像素点是否可能是极值
//调用函数pixval32f获取图像dog_pyr[o][i]的第r行第c列的点的坐标值，
//然后调用ABS宏求其绝对值
if( ABS( pixval32f( dog_pyr[o][i], r, c ) ) > prelim_contr_thr )
//通过在尺度空间中将一个像素点的值与其周围3*3*3邻域内的点比较来
//决定此点是否极值点(极大值或极小都行)
if( is_extremum( dog_pyr, o, i, r, c ) )//若是极值点
{
//由于极值点的检测是在离散空间中进行的，所以检测到的极值点并
//不一定是真正意义上的极值点
//因为真正的极值点可能位于两个像素之间，而在离散空间中只能精
//确到坐标点精度上
//通过亚像素级插值进行极值点精确定位(修正极值点坐标)，并去除
//低对比度的极值点，将修正后的特征点组成feature结构返回
feat = interp_extremum(dog_pyr, o, i, r, c, intvls, contr_thr);
//返回值非空，表明此点已被成功修正
if( feat )
{
//调用宏feat_detection_data来提取参数feat中的feature_data成员
//并转换为detection_data类型的指针
ddata = feat_detection_data( feat );
//去除边缘响应，即通过计算主曲率比值判断某点是否边缘点,
//返回值为0表示不是边缘点，可做特征点
if( ! is_too_edge_like( dog_pyr[ddata->octv][ddata->intvl],
ddata->r, ddata->c, curv_thr ) )
{
cvSeqPush( features, feat );//向特征点序列features末尾插
//入新检测到的特征点feat
}
else
free( ddata );
free( feat );
}
}
return features;//返回特征点序列
}
（5.2.1）CvSeq的理解：
说明：动态结构序列CvSeq是所有OpenCv动态数据结构的基础。有两种类型的序列：稠密序列，稀疏序列：
其中一种是——稠密序列都派生自CvSeq，他们用来代表可扩展的一维数组 — 向量、栈、队列和双端队列。数据间不存在空隙（连续存储）。如果元素元素从序列中间被删除或插入新的元素到序列，那么此元素后边的相关元素全部被移动。
total表示稠密序列的元素个数，或者稀疏序列被分配的节点数。elem_size表示序列中每个元素占用的字节数。block_max是最近一个内存的最大边界指针。ptr表示当写指针。delta_elems表示序列间隔尺寸。storage指向序列存储的内存块的指针。free_blocks表示空的块列表。first指向第一个序列块。
参看博客：1）OpenCV——CvSeq动态结构序列：https://www.cnblogs.com/farewell-farewell/p/5999908.html
2）CvSeq的理解：https://blog.csdn.net/weijianmeng/article/details/7173560
（5.2.2）cvCreateSeq的理解：
说明：CvSeq *cvCreateSeq(int seq_flags, int header_size, int elem_size, CvMemStorage *storage);
header_size：序列头的大小，通常为sizeof(CvSeq)。
elem_size：存储元素的大小。
storage：内存存储器，添加元素的时候，就会从内存存储器申请空间。
参看博客：opencv创建序列cvCreateSeq与插入元素cvSeqPush的运用：https://www.cnblogs.com/farewell-farewell/p/5999908.html
（5.2.3）for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)//第r行
for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)//第c列
中为什么需要设置SIFT_IMG_BORDER指明边界宽度：
答：只检测边界线以内的极值点，因为边界容易受到噪声的干扰而不稳定。
（5.2.4）pixval32f(dog_pyr[o][i])函数功能？ABS(pixval32f(dog_pyr[o][i]))为什么加绝对值？为什么需要判断高斯差分金字塔（DOG）中某个像素小于0的情况？
答：utils.c中定义了这个内部函数：
static inline float pixval32f( IplImage* img, int r, int c )
{
return ( (float*)(img->imageData + img->widthStep*r) )[c];
}
具体怎么访问图像元素？
参看博客：1） opencv——访问图像元素(imagedata widthstep)：https://blog.csdn.net/xiaofeilong321/article/details/13287697
2）opencv学习笔记（八）IplImage* 访问图像像素的值：https://www.cnblogs.com/codingmengmeng/p/6559724.html
（5.2.5）is_extremum代码及说明：
/*通过在尺度空间中将一个像素点的值与其周围3*3*3邻域内的点比较来决定此点是否极值点
(极大值或极小都行)
参数：
dog_pyr：高斯差分金字塔
octv：像素点所在的组
intvl：像素点所在的层
r：像素点所在的行
c：像素点所在的列
返回值：若指定的像素点是极值点(极大值或极小值)，返回1；否则返回0*/
static int is_extremum( IplImage*** dog_pyr, int octv, int intvl, int r, int c )
{
//调用函数pixval32f获取图像dog_pyr[octv][intvl]的第r行第c列的点的坐标值
float val = pixval32f( dog_pyr[octv][intvl], r, c );
int i, j, k;
//检查是否最大值
if( val > 0 )
{
for( i = -1; i <= 1; i++ )//层
for( j = -1; j <= 1; j++ )//行
for( k = -1; k <= 1; k++ )//列
if( val < pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )
return 0; //不是极值点，退出此函数
}
//检查是否最小值
else
{
for( i = -1; i <= 1; i++ )//层
for( j = -1; j <= 1; j++ )//行
for( k = -1; k <= 1; k++ )//列
if( val > pixval32f( dog_pyr[octv][intvl+i], r + j, c + k ) )
return 0;
}
return 1; //是极值点，返回1
}
（5.2.6）interp_extremum代码及说明：
剔除不稳定点，精确定位关键点位置
/*通过亚像素级插值进行极值点精确定位(修正极值点坐标)，并去除低对比度的极值点，
将修正后的特征点组成feature结构返回
参数：
dog_pyr：高斯差分金字塔
octv：像素点所在的组
intvl：像素点所在的层
r：像素点所在的行
c：像素点所在的列
intvls：每组的层数
contr_thr：对比度阈值，针对归一化后的图像，用来去除不稳定特征
返回值：返回经插值修正后的特征点(feature类型)；若经有限次插值依然无法精确到理想情况
或者该点对比度过低，返回NULL*/
static struct feature* interp_extremum( IplImage*** dog_pyr, int octv, int intvl,
int r, int c, int intvls, double contr_thr ) //intvls=3; contr_thr=0.04
{
struct feature* feat;//修正后的特征点
struct detection_data* ddata;//与特征检测有关的结构，存在feature结构的feature_data成员中
double xi, xr, xc, contr;//xi,xr,xc分别为亚像素的intvl(层),row(y),col(x)方向上的
//增量(偏移量)
int i = 0;//插值次数
//SIFT_MAX_INTERP_STEPS指定了关键点的最大插值次数，即最多修正多少次，默认是5
while( i < SIFT_MAX_INTERP_STEPS )
{
//进行一次极值点差值，计算σ(层方向,intvl方向)，y，x方向上的子像素偏移量(增量)
interp_step( dog_pyr, octv, intvl, r, c, &xi, &xr, &xc );
//若在任意方向上的偏移量大于0.5时，意味着差值中心已经偏移到它的临近点上，
//所以必须改变当前关键点的位置坐标
if( ABS( xi ) < 0.5 && ABS( xr ) < 0.5 && ABS( xc ) < 0.5 )//若三方向上偏移量
//都小于0.5，表示已经够精确，则不用继续插值
break;
//修正关键点的坐标，x,y,σ三方向上的原坐标加上偏移量取整(四舍五入)
c += cvRound( xc );//x坐标修正
r += cvRound( xr );//y坐标修正
intvl += cvRound( xi );//σ方向，即层方向
//若坐标修正后超出范围，则结束插值，返回NULL
if( intvl < 1 || //层坐标插之后越界
intvl > intvls ||
c < SIFT_IMG_BORDER || //行列坐标插之后到边界线内
r < SIFT_IMG_BORDER ||
c >= dog_pyr[octv][0]->width - SIFT_IMG_BORDER ||
r >= dog_pyr[octv][0]->height - SIFT_IMG_BORDER )
{
return NULL;
}
i++;
}
//若经过SIFT_MAX_INTERP_STEPS次插值后还没有修正到理想的精确位置，则返回NULL，
//即舍弃此极值点
if( i >= SIFT_MAX_INTERP_STEPS )
return NULL;
//计算被插值点的对比度：D + 0.5 * dD^T * X
contr = interp_contr( dog_pyr, octv, intvl, r, c, xi, xr, xc );
//此处与书上的不太一样，RobHess此处的阈值使用的是0.04/S
if( ABS( contr ) < contr_thr / intvls )//若该点对比度过小，舍弃，返回NULL
return NULL;
//为一个特征点feature结构分配空间并初始化，返回特征点指针
feat = new_feature();
//调用宏feat_detection_data来提取参数feat中的feature_data成员并转换为
//detection_data类型的指针
ddata = feat_detection_data( feat );
//将修正后的坐标赋值给特征点feat
//原图中特征点的x坐标，因为第octv组中的图的尺寸比原图小2^octv倍，
//所以坐标值要乘以2^octv,最后针对的是原图的匹配
feat->img_pt.x = feat->x = ( c + xc ) * pow( 2.0, octv );
//原图中特征点的y坐标，因为第octv组中的图的尺寸比原图小2^octv倍，
//所以坐标值要乘以2^octv
feat->img_pt.y = feat->y = ( r + xr ) * pow( 2.0, octv );
ddata->r = r;//特征点所在的行
ddata->c = c;//特征点所在的列
ddata->octv = octv;//高斯差分金字塔中，特征点所在的组
ddata->intvl = intvl;//高斯差分金字塔中，特征点所在的组中的层
ddata->subintvl = xi;//特征点在层方向(σ方向,intvl方向)上的亚像素偏移量
return feat;//返回特征点指针
}
（5.2.6.1）cvRound()说明：
答：函数cvRound，cvFloor，cvCeil 都是用一种舍入的方法将输入浮点数转换成整数：
cvRound()：返回跟参数最接近的整数值，即四舍五入；
cvFloor()：返回不大于参数的最大整数值，即向下取整；
cvCeil()：返回不小于参数的最小整数值，即向上取整；
参看:【杂谈opencv】OpenCV中的cvRound()、cvFloor()、 cvCeil()函数讲解——https://blog.csdn.net/sinat_36264666/article/details/78849125
（5.2.6.2）interp_step代码及说明：
答：
/*进行一次极值点差值，计算x，y，σ方向(层方向)上的子像素偏移量(增量)
参数：
dog_pyr：高斯差分金字塔
octv：像素点所在的组
intvl：像素点所在的层
r：像素点所在的行
c：像素点所在的列
xi：输出参数，层方向上的子像素增量(偏移)
xr：输出参数，y方向上的子像素增量(偏移)
xc：输出参数，x方向上的子像素增量(偏移)*/
static void interp_step( IplImage*** dog_pyr, int octv, int intvl, int r, int c,
double* xi, double* xr, double* xc )
{
CvMat* dD, * H, * H_inv, X;
double x[3] = { 0 };
//在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数，结果存放在列向量dD中
dD = deriv_3D( dog_pyr, octv, intvl, r, c );
//在DoG金字塔中计算某点的3*3海森矩阵
H = hessian_3D( dog_pyr, octv, intvl, r, c );
H_inv = cvCreateMat( 3, 3, CV_64FC1 );//海森矩阵的逆阵
cvInvert( H, H_inv, CV_SVD );
cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );
//X = - H^(-1) * dD，H的三个元素分别是x,y,σ方向上的偏移量(具体见SIFT算法说明)
cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );
cvReleaseMat( &dD );
cvReleaseMat( &H );
cvReleaseMat( &H_inv );
*xi = x[2];//σ方向(层方向)偏移量
*xr = x[1];//y方向上偏移量
*xc = x[0];//x方向上偏移量
}
（5.2.6.2.1）deriv_3D代码及说明：
答：
/*在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数
参数：
dog_pyr：高斯差分金字塔
octv：像素点所在的组
intvl：像素点所在的层
r：像素点所在的行
c：像素点所在的列
返回值：返回3个偏导数组成的列向量{ dI/dx, dI/dy, dI/ds }^T*/
static CvMat* deriv_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c )
{
CvMat* dI;
double dx, dy, ds;
//求差分来代替偏导，这里是用的隔行求差取中值的梯度计算方法
//求x方向上的差分来近似代替偏导数
dx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) -
pixval32f( dog_pyr[octv][intvl], r, c-1 ) ) / 2.0;
//求y方向上的差分来近似代替偏导数
dy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) -
pixval32f( dog_pyr[octv][intvl], r-1, c ) ) / 2.0;
//求层间的差分来近似代替尺度方向上的偏导数
ds = ( pixval32f( dog_pyr[octv][intvl+1], r, c )
pixval32f( dog_pyr[octv][intvl-1], r, c ) ) / 2.0;
//组成列向量
dI = cvCreateMat( 3, 1, CV_64FC1 );
cvmSet( dI, 0, 0, dx );
cvmSet( dI, 1, 0, dy );
cvmSet( dI, 2, 0, ds );
return dI;
}
（5.2.6.2.1.1）问题issue：差分代替偏导？一阶二阶都使用吗？怎么操作？
差分代替偏导！一阶情况：
一阶差分代替一阶偏导，即梯度，[-1，0，1]
行方向[-1，0，1]
-1
列方向[ 0 ]
1
二阶情况：
参看：图像处理中的一阶偏导数和二阶偏导数——https://www.docin.com/p-749102193.html
（5.2.6.2.1.2）cvmSet说明
答：opencv中cvmSet为逐点赋值
cvmSet( dI, 0, 0, dx ); dx
cvmSet( dI, 1, 0, dy ); dI = [ dy ]
cvmSet( dI, 2, 0, ds ); ds
（5.2.6.2.2）hessian_3D代码及说明：
答：//差分代替偏导，二阶情况
/*在DoG金字塔中计算某点的3*3海森矩阵
/ Ixx Ixy Ixs \
| Ixy Iyy Iys |
\ Ixs Iys Iss /
参数：
dog_pyr：高斯差分金字塔
octv：像素点所在的组
intvl：像素点所在的层
r：像素点所在的行
c：像素点所在的列
返回值：返回3*3的海森矩阵
*/
static CvMat* hessian_3D( IplImage*** dog_pyr, int octv, int intvl, int r, int c )
{
CvMat* H;
double v, dxx, dyy, dss, dxy, dxs, dys;
v = pixval32f( dog_pyr[octv][intvl], r, c );//该点的像素值
//用差分近似代替倒数(具体公式见各种梯度的求法)
//dxx = f(i+1,j) - 2f(i,j) + f(i-1,j)
//dyy = f(i,j+1) - 2f(i,j) + f(i,j-1)
dxx = ( pixval32f( dog_pyr[octv][intvl], r, c+1 ) +
pixval32f( dog_pyr[octv][intvl], r, c-1 ) - 2 * v );
dyy = ( pixval32f( dog_pyr[octv][intvl], r+1, c ) +
pixval32f( dog_pyr[octv][intvl], r-1, c ) - 2 * v );
dss = ( pixval32f( dog_pyr[octv][intvl+1], r, c ) +
pixval32f( dog_pyr[octv][intvl-1], r, c ) - 2 * v );
dxy = ( pixval32f( dog_pyr[octv][intvl], r+1, c+1 ) -
pixval32f( dog_pyr[octv][intvl], r+1, c-1 ) -
pixval32f( dog_pyr[octv][intvl], r-1, c+1 ) +
pixval32f( dog_pyr[octv][intvl], r-1, c-1 ) ) / 4.0;
dxs = ( pixval32f( dog_pyr[octv][intvl+1], r, c+1 ) -
pixval32f( dog_pyr[octv][intvl+1], r, c-1 ) -
pixval32f( dog_pyr[octv][intvl-1], r, c+1 ) +
pixval32f( dog_pyr[octv][intvl-1], r, c-1 ) ) / 4.0;
dys = ( pixval32f( dog_pyr[octv][intvl+1], r+1, c ) -
pixval32f( dog_pyr[octv][intvl+1], r-1, c ) -
pixval32f( dog_pyr[octv][intvl-1], r+1, c ) +
pixval32f( dog_pyr[octv][intvl-1], r-1, c ) ) / 4.0;
//组成海森矩阵
H = cvCreateMat( 3, 3, CV_64FC1 );
cvmSet( H, 0, 0, dxx );
cvmSet( H, 0, 1, dxy );
cvmSet( H, 0, 2, dxs );
cvmSet( H, 1, 0, dxy );
cvmSet( H, 1, 1, dyy );
cvmSet( H, 1, 2, dys );
cvmSet( H, 2, 0, dxs );
cvmSet( H, 2, 1, dys );
cvmSet( H, 2, 2, dss );
return H;
}
（5.2.6.2.3）cvInvert函数说明：
答：
double cvInvert(//矩阵取逆
const CvArr* src,//目标矩阵
CvArr* dst,//结果矩阵
int method = CV_LU//逆运算方法
);
其中method有
方法的参数值 含义
CV_LU 高斯消去法
CV_SVD 奇异值分解
CV_SVD_SYM 对称矩阵的SVD
参看：《学习opencv》笔记——矩阵和图像操作——cvInRange,cvInRangeS,cvInvert and cvMahalonobis——https://blog.csdn.net/zhurui_idea/article/details/28630589
（5.2.6.2.4）问题issue：cvInitMatHeader的一些问题：
答：
1. cvInitMatHeader的用法很古怪。第一个参数必须是CvMat格式的，官方的文档是CvMat *mat，但是这里要注意，mat必须是初始化过的，单独定义一个指针如CvMat *data1；把data1带入cvInitMatHeader函数后，编译器会报错，显示没有初始化！下面的代码才是正确的
2.下面的代码我在使用的时候想生成两个mat，然后被另一个函数f调用，发现f调用时失败的，用下面代码生成的数据被释放掉了。所以这种用法只能在一个函数体内使用，跨函数是不行的。
CvMat data1;
CvMat responses1;
//将数组中数据转化成opencv支持的mat格式
cvInitMatHeader(&data1, total, var_count, CV_32FC1, pData);
cvInitMatHeader(&responses1, total, 1, CV_32SC1, pRespon);
参看：1）【OpenCV学习】矩阵cvInitMatHeader和cvCreateMat：https://blog.csdn.net/cc1949/article/details/22476251
2）cvInitMatHeader的一些问题：https://www.xuebuyuan.com/1050480.html
3）【OpenCV学习】矩阵cvInitMatHeader和cvCreateMat：举例说明看——https://blog.csdn.net/cc1949/article/details/22476251
（5.2.6.2.5）cvGEMM函数说明：
答：void cvGEMM( const CvArr* src1, const CvArr* src2, double alpha, const CvArr* src3, double beta, CvArr* dst, int tABC=0 );
这是通用矩阵乘法，其中各个参数表示：
src1：第一输入数组
src2：第二输入数组
alpha：系数
src3“第三输入数组（偏移量），如果没有偏移量，可以为空（NULL）
beta：表示偏移量的系数
dst：输出数组
tABC：转置操作标志，可以是0。当为0时，没有转置。
参看博客：cvGEMM、cvMatMul和cvMatMulAdd的定义——https://blog.csdn.net/liulianfanjianshi/article/details/11737921
（5.2.6.3）interp_contr代码及说明：
答：
/*计算被插值点的对比度：D + 0.5 * dD^T * X
参数：
dog_pyr：高斯差分金字塔
octv：像素点所在的组
intvl：像素点所在的层
r：像素点所在的行
c：像素点所在的列
xi：层方向上的子像素增量
xr：y方向上的子像素增量
xc：x方向上的子像素增量
返回值：插值点的对比度*/
static double interp_contr( IplImage*** dog_pyr, int octv, int intvl, int r,
int c, double xi, double xr, double xc )
{
CvMat* dD, X, T;
double t[1], x[3] = { xc, xr, xi };
//偏移量组成的列向量X，其中是x，y，σ三方向上的偏移量
cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );
//矩阵乘法的结果T，是一个数值
cvInitMatHeader( &T, 1, 1, CV_64FC1, t, CV_AUTOSTEP );
//在DoG金字塔中计算某点的x方向、y方向以及尺度方向上的偏导数，结果存放在列向量dD中
dD = deriv_3D( dog_pyr, octv, intvl, r, c );
//矩阵乘法：T = dD^T * X
cvGEMM( dD, &X, 1, NULL, 0, &T, CV_GEMM_A_T );
cvReleaseMat( &dD );
//返回计算出的对比度值：D + 0.5 * dD^T * X (具体公式推导见SIFT算法说明)
return pixval32f( dog_pyr[octv][intvl], r, c ) + t[0] * 0.5;
}
（5.2.6.4）new_feature()代码及说明：
答：
/*为一个feature结构分配空间并初始化
返回值：初始化完成的feature结构的指针*/
static struct feature* new_feature( void )
{
struct feature* feat;//特征点指针
struct detection_data* ddata;//与特征检测相关的结构
feat = malloc( sizeof( struct feature ) );//分配空间
memset( feat, 0, sizeof( struct feature ) );//清零
ddata = malloc( sizeof( struct detection_data ) );
memset( ddata, 0, sizeof( struct detection_data ) );
feat->feature_data = ddata;//将特征检测相关的结构指针赋值给特征点的feature_data成员
feat->type = FEATURE_LOWE;//默认是LOWE类型的特征点
return feat;
}
（5.2.6.5）pow()说明：
答：C 库函数 double pow(double x, double y) 返回 x 的 y 次幂，即 x^y。
参看:C 库函数 - pow()——https://www.runoob.com/cprogramming/c-function-pow.html
（5.2.7）is_too_edge_like代码及说明：
/*去除边缘响应，即通过计算主曲率比值判断某点是否边缘点
参数：
dog_img：此特征点所在的DoG图像
r：特征点所在的行
c：特征点所在的列
cur_thr：主曲率比值的阈值，用来去除边缘特征
返回值：0：此点是非边缘点；1：此点是边缘点*/
static int is_too_edge_like( IplImage* dog_img, int r, int c, int curv_thr )
{
double d, dxx, dyy, dxy, tr, det;
/*某点的主曲率与其海森矩阵的特征值成正比，为了避免直接计算特征值，这里只考虑
特征值的比值可通过计算海森矩阵的迹tr(H)和行列式det(H)来计算特征值的比值
设a是海森矩阵的较大特征值，b是较小的特征值，有a = r*b，r是大小特征值的比值
tr(H) = a + b; det(H) = a*b;
tr(H)^2 / det(H) = (a+b)^2 / ab = (r+1)^2/r
r越大，越可能是边缘点；伴随r的增大，(r+1)^2/r 的值也增大，所以可通过(r+1)^2/r 判断
主曲率比值是否满足条件*/
/* principal curvatures are computed using the trace and det of Hessian */
d = pixval32f(dog_img, r, c);//调用函数pixval32f获取图像dog_img的第r行第c列的点的坐标值
//用差分近似代替偏导，求出海森矩阵的几个元素值
/* / dxx dxy \
\ dxy dyy / */
dxx = pixval32f( dog_img, r, c+1 ) + pixval32f( dog_img, r, c-1 ) - 2 * d;
dyy = pixval32f( dog_img, r+1, c ) + pixval32f( dog_img, r-1, c ) - 2 * d;
dxy = ( pixval32f(dog_img, r+1, c+1) - pixval32f(dog_img, r+1, c-1) -
pixval32f(dog_img, r-1, c+1) + pixval32f(dog_img, r-1, c-1) ) / 4.0;
tr = dxx + dyy;//海森矩阵的迹
det = dxx * dyy - dxy * dxy;//海森矩阵的行列式
//若行列式为负，表明曲率有不同的符号，去除此点
/* negative determinant -> curvatures have different signs; reject feature */
if( det <= 0 )
return 1;//返回1表明此点是边缘点
//通过式子：(r+1)^2/r 判断主曲率的比值是否满足条件，若小于阈值，表明不是边缘点
if( tr * tr / det < ( curv_thr + 1.0 )*( curv_thr + 1.0 ) / curv_thr )
return 0;//不是边缘点
return 1;//是边缘点
}
（5.2.8）cvSeqPush代码及说明：
函数原型： char * cvSeqPush(CvSeq *seq, void *element = NULL);//添加元素到序列的尾部
参看：opencv创建序列cvCreateSeq与插入元素cvSeqPush的运用：https://blog.csdn.net/gdut2015go/article/details/46494677
（5.3）calc_feature_scales代码及说明：
/*计算特征点序列中每个特征点的尺度
参数：
features：特征点序列
sigma：初始高斯平滑参数，即初始尺度
intvls：尺度空间中每组的层数*/
static void calc_feature_scales( CvSeq* features, double sigma, int intvls )
{
struct feature* feat;
struct detection_data* ddata;
double intvl;
int i, n;
n = features->total;//总的特征点个数
//遍历特征点
for( i = 0; i < n; i++ )
{
//调用宏，获取序列features中的第i个元素，并强制转换为struct feature类型
feat = CV_GET_SEQ_ELEM( struct feature, features, i );
//调用宏feat_detection_data来提取参数feat中的feature_data成员并转换为
//detection_data类型的指针
ddata = feat_detection_data( feat );
//特征点所在的层数ddata->intvl加上特征点在层方向上的亚像素偏移量，得到
//特征点的较为精确的层数
intvl = ddata->intvl + ddata->subintvl;
//计算特征点的尺度(公式见SIFT算法说明)，并赋值给scl成员
feat->scl = sigma * pow( 2.0, ddata->octv + intvl / intvls );
//计算特征点所在的组的尺度，给detection_data的scl_octv成员赋值
ddata->scl_octv = sigma * pow( 2.0, intvl / intvls );
}
}
（5.3.1）CV_GET_SEQ_ELEM函数说明：
答案：用法：从所给序列中取出元素的地址，注意：得到的是地址，即指针
所以，关键在于序列中存放的是那种类型的数据，若存放的为地址，那用这个宏得到的就是指针的指针。
例程：
1.序列中存放的为CvPoint
CvPoint pt=*CV_GET_SEQ_ELEM(CvPoint,hull,i);
2.序列中存放的为CvPoint*，即指针
CvPoint* pt=*CV_GET_SEQ_ELEM(CvPoint*,hull,i);
或者
CvPoint pt=**CV_GET_SEQ_ELEM(CvPoint*,hull,i);
参看：关于宏CV_GET_SEQ_ELEM——https://blog.csdn.net/u013089125/article/details/20126215
（5.4）adjust_for_img_dbl代码及说明：
/*将特征点序列中每个特征点的坐标减半(当设置了将图像放大为原图的2倍时，特征点检测
完之后调用)
参数：
features：特征点序列*/
static void adjust_for_img_dbl( CvSeq* features )
{
struct feature* feat;
int i, n;
n = features->total;//总的特征点个数
//遍历特征点
for( i = 0; i < n; i++ )
{
//调用宏，获取序列features中的第i个元素，并强制转换为struct feature类型
feat = CV_GET_SEQ_ELEM( struct feature, features, i );
//将特征点的x,y坐标和尺度都减半
feat->x /= 2.0;
feat->y /= 2.0;
feat->scl /= 2.0;
feat->img_pt.x /= 2.0;
feat->img_pt.y /= 2.0;
}
}
附：
1.问题描述：C版本的cvGetTickFrequency()函数和C++版本的getTickFrequency()计算时间的公式不一样?
解决方案：参看博客：https://blog.csdn.net/chaipp0607/article/details/71056580
注意：OpenCV C版本的cvGetTickFrequency()函数和C++版本的getTickFrequency()的单位不一样
posted on 2019-06-03 17:29 Alliswell_WP 阅读(1888) 评论(0) 编辑 收藏 举报
