关于SIFT的原理已经有很多大牛的博客上做了解析，本文重点将以Rob
Hess等人用C实现的代码做解析，结合代码SIFT原理会更容易理解。一些难理解点的用了☆标注。
欢迎大家批评指正！
SIFT（Scale-invariant feature transform）即尺度不变特征转换，提取的局部特征点具有尺度不变性，且对于旋转，亮度，噪声等有很高的稳定性。
SIFT特征点提取
本文将以下函数为参照顺序介绍SIFT特征点提取与描述方法。
1.图像预处理
2.构建高斯金字塔（不同尺度下的图像）
3.生成DOG尺度空间
4.关键点搜索与定位
5.计算特征点对应原图的位置
6.为特征点分配方向角
7.构建特征描述子
/**
Finds SIFT features in an image using user-specified parameter values. All
detected features are stored in the array pointed to by \a feat.
*/
int _sift_features( IplImage* img, struct feature** feat, int intvls,
double sigma, double contr_thr, int curv_thr,
int img_dbl, int descr_width, int descr_hist_bins )
{
IplImage* init_img;
IplImage*** gauss_pyr, *** dog_pyr;
CvMemStorage* storage;
CvSeq* features;
int octvs, i, n = 0;
/* check arguments */
if( ! img )
fatal_error( "NULL pointer error, %s, line %d", __FILE__, __LINE__ );
if( ! feat )
fatal_error( "NULL pointer error, %s, line %d", __FILE__, __LINE__ );
/* build scale space pyramid; smallest dimension of top level is ~4 pixels */
init_img = create_init_img( img, img_dbl, sigma ); //对进行图片预处理
octvs = log( MIN( init_img->width, init_img->height ) ) / log(2) - 2; //计算高斯金字塔的组数(octave),同时保证顶层至少有4个像素点
gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma ); //建立高斯金字塔
dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls ); //DOG尺度空间
storage = cvCreateMemStorage( 0 );
features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr, //极值点检测，并去除不稳定特征点
curv_thr, storage );
calc_feature_scales( features, sigma, intvls ); //计算特征点所在的尺度
if( img_dbl )
adjust_for_img_dbl( features ); //如果图像初始被扩大了2倍，所有坐标与尺度要除以2
calc_feature_oris( features, gauss_pyr ); //计算特征点的方向角
compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins );//计算特征描述子
/* sort features by decreasing scale and move from CvSeq to array */
cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL ); //对特征点按尺度排序
n = features->total;
*feat = calloc( n, sizeof(struct feature) );
*feat = cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ );
for( i = 0; i < n; i++ )
{
free( (*feat)[i].feature_data );
(*feat)[i].feature_data = NULL;
}
cvReleaseMemStorage( &storage );
cvReleaseImage( &init_img );
release_pyr( &gauss_pyr, octvs, intvls + 3 );
release_pyr( &dog_pyr, octvs, intvls + 2 );
return n;
}
—————————————————————————————————————————————————————
1.图像预处理
/************************ Functions prototyped here **************************/
/*
Converts an image to 8-bit grayscale and Gaussian-smooths it. The image is
optionally doubled in size prior to smoothing.
@param img input image
@param img_dbl if true, image is doubled in size prior to smoothing
@param sigma total std of Gaussian smoothing
*/
static IplImage* create_init_img( IplImage* img, int img_dbl, double sigma )
{
IplImage* gray, * dbl;
double sig_diff;
gray = convert_to_gray32( img ); //转换为32位灰度图
if( img_dbl ) // 图像被放大二倍
{
sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4 ); // sigma = 1.6 , SIFT_INIT_SIGMA = 0.5 lowe认为图像在尺度0.5下最清晰
dbl = cvCreateImage( cvSize( img->width*2, img->height*2 ),
IPL_DEPTH_32F, 1 );
cvResize( gray, dbl, CV_INTER_CUBIC ); //双三次插值方法 放大图像
cvSmooth( dbl, dbl, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff ); //高斯平滑
cvReleaseImage( &gray );
return dbl;
}
else
{
sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA );
cvSmooth( gray, gray, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff ); // 高斯平滑
return gray;
}
}
lowe建议把初始图像放大二倍，可以得到更多的特征点，提取到更多细节，并且认为图像在尺度σ = 0.5时图像最清晰，初始高斯尺度为σ = 1.6。
☆第19行因为图像被放大二倍，此时σ = 1.0 。因为对二倍化后的图像平滑是在σ = 0.5 上叠加的高斯模糊，
所以有模糊系数有sig_diff = sqrt (sigma *sigma - 0.5*0.5*4)=sqrt(1.6*1.6 -1) ；
2.构建高斯金字塔
构建高斯金字塔过程即构建出图像在不同尺度上图像，提取到的特征点可有具有尺度不变性。
图像的尺度空间L(x,y,σ)可以用一个高斯函数G(x,y,σ)与图像I(x,y)卷积产生，即L(x,y,σ) = G(x,y,σ) * I(x,y)
其中二维高斯核的计算为 
。
☆不同的尺度空间即用不同的高斯核函数平滑图像， 平滑系数越大，图像越模糊。即模拟出动物的视觉效果，因为事先不知道物体的大小，在不同的尺度下，图像的细节会表现的不同。当尺度由小变大的过程中，是一个细节逐步简化的过程，图像中特征不够明显的物体，就模糊的多了，而有些物体还可以看得到大致的轮廓。所以要在不同尺度下，观察物体的尺度响应，提取到的特征才能具有尺度不变性。
SIFT算法采用高斯金字塔实现连续的尺度空间的图像。金字塔共分为O(octave)组，每组有S(intervals)层 ，下一组是由上一组隔点采样得到（即降2倍分辨率），这是为了减轻卷积运算的工作量。
构建高斯金字塔（octave = 5, intervals +3=6）：
全部空间尺度为：
☆1.这个尺度因子都是在原图上进行的，而在算法实现过程中，采用高斯平滑是在上一层图像上再叠加高斯平滑，即我们在程序中看到的平滑因子为
Eg. 在第一层上为了得到kσ的高斯模糊图像，可以在原图上直接采用kσ平滑，也可以在上一层图像上（已被σ平滑）的图像上采用平滑因子为
平滑图像，效果是一样的。
☆2.我们在源码上同时也没有看到组间的2倍的关系，实际在对每组的平滑因子都是一样的，2倍的关系是由于在降采样的过程中产生的，第二层的第一张图是由第一层的平滑因子为2σ的图像上（即倒数第三张）降采样得到，此时图像平滑因子为σ，所以继续采用以上的平滑因子。但在原图上看，形成了全部的空间尺度。
☆3.每组（octave）有S+3层图像，是由于在DOG尺度空间上寻找极值点的方法是在一个立方体内进行，即上下层比较，所以不在DOG空间的第一层与最后一层寻找，即DOG需要S+2层图像，由于DOG尺度空间是由高斯金字塔相邻图像相减得到，即每组需要S+3层图像。
/*
Builds Gaussian scale space pyramid from an image
@param base base image of the pyramid
@param octvs number of octaves of scale space
@param intvls number of intervals per octave
@param sigma amount of Gaussian smoothing per octave
@return Returns a Gaussian scale space pyramid as an octvs x (intvls + 3)
array
给定组数(octave)和层数(intvls)，以及初始平滑系数sigma,构建高斯金字塔
返回的每组中层数为intvls+3
*/
static IplImage*** build_gauss_pyr( IplImage* base, int octvs,
int intvls, double sigma )
{
IplImage*** gauss_pyr;
const int _intvls = intvls; // lowe 采用了每组层数(intvls)为 3
// double sig_total, sig_prev;
double k;
int i, o;
double *sig = (double *)malloc(sizeof(int)*(_intvls+3)); //存储每组的高斯平滑因子，每组对应的平滑因子都相同
gauss_pyr = calloc( octvs, sizeof( IplImage** ) );
for( i = 0; i < octvs; i++ )
gauss_pyr[i] = calloc( intvls + 3, sizeof( IplImage *) );
/*
precompute Gaussian sigmas using the following formula:
\sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2
sig[i] is the incremental sigma value needed to compute
the actual sigma of level i. Keeping track of incremental
sigmas vs. total sigmas keeps the gaussian kernel small.
*/
k = pow( 2.0, 1.0 / intvls ); // k = 2^(1/S)
sig[0] = sigma;
sig[1] = sigma * sqrt( k*k- 1 );
for (i = 2; i < intvls + 3; i++)
sig[i] = sig[i-1] * k; //每组对应的平滑因子为 σ , sqrt(k^2 -1)* σ, sqrt(k^2 -1)* kσ , ...
for( o = 0; o < octvs; o++ )
for( i = 0; i < intvls + 3; i++ )
{
if( o == 0 && i == 0 )
gauss_pyr[o][i] = cvCloneImage(base); //第一组，第一层为原图
/* base of new octvave is halved image from end of previous octave */
else if( i == 0 )
gauss_pyr[o][i] = downsample( gauss_pyr[o-1][intvls] ); //第一层图像由上一层倒数第三张隔点采样得到
/* blur the current octave's last image to create the next one */
else
{
gauss_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i-1]),
IPL_DEPTH_32F, 1 );
cvSmooth( gauss_pyr[o][i-1], gauss_pyr[o][i],
CV_GAUSSIAN, 0, 0, sig[i], sig[i] ); //高斯平滑
}
}
return gauss_pyr;
}
3.生成DOG尺度空间
Lindeberg发现高斯差分函数（Difference of Gaussian ，简称DOG算子）与尺度归一化的高斯拉普拉斯函数
非常近似，且
差分近似：
lowe建议采用相邻尺度的图像相减来获得高斯差分图像D(x,y,σ)来近似LOG来进行极值检测。
D(x,y,σ) = G(x,y,kσ)*I(x,y)-G(x,y,σ)*I(x,y)
=L(x,y,kσ) - L(x,y,σ)
对高斯金字塔的每组内相邻图像相减，形成DOG尺度空间，这时DOG中每组有S+2层图像
static IplImage*** build_dog_pyr( IplImage*** gauss_pyr, int octvs, int intvls )
{
IplImage*** dog_pyr;
int i, o;
dog_pyr = calloc( octvs, sizeof( IplImage** ) );
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_pyr[o][i] = cvCreateImage( cvGetSize(gauss_pyr[o][i]),
IPL_DEPTH_32F, 1 );
cvSub( gauss_pyr[o][i+1], gauss_pyr[o][i], dog_pyr[o][i], NULL ); //相邻两层图像相减，结果放在dog_pyr数组内
}
return dog_pyr;
}
4.关键点搜索与定位
在DOG尺度空间上，首先寻找极值点，插值处理，找到准确的极值点坐标，再排除不稳定的特征点（边界点）
/*
Detects features at extrema in DoG scale space. Bad features are discarded
based on contrast and ratio of principal curvatures.
@return Returns an array of detected features whose scales, orientations,
and descriptors are yet to be determined.
*/
static CvSeq* scale_space_extrema( IplImage*** dog_pyr, int octvs, int intvls,
double contr_thr, int curv_thr,
CvMemStorage* storage )
{
CvSeq* features;
double prelim_contr_thr = 0.5 * contr_thr / intvls; //极值比较前的阈值处理
struct feature* feat;
struct detection_data* ddata;
int o, i, r, c;
features = cvCreateSeq( 0, sizeof(CvSeq), sizeof(struct feature), storage );
for( o = 0; o < octvs; o++ ) //对DOG尺度空间上，遍历从第二层图像开始到倒数第二层图像上，每个像素点
for( i = 1; i <= intvls; i++ )
for(r = SIFT_IMG_BORDER; r < dog_pyr[o][0]->height-SIFT_IMG_BORDER; r++)
for(c = SIFT_IMG_BORDER; c < dog_pyr[o][0]->width-SIFT_IMG_BORDER; c++)
/* perform preliminary check on contrast */
if( ABS( pixval32f( dog_pyr[o][i], r, c ) ) > prelim_contr_thr ) // 排除像素值小于阈值prelim_contr_thr的点，提高稳定性
if( is_extremum( dog_pyr, o, i, r, c ) ) //与周围26个像素值比较，是否极大值或者极小值点
{
feat = interp_extremum(dog_pyr, o, i, r, c, intvls, contr_thr); //插值处理，找到准确的特征点坐标
if( feat )
{
ddata = feat_detection_data( feat );
if( ! is_too_edge_like( dog_pyr[ddata->octv][ddata->intvl], //根据Hessian矩阵 判断是否为边缘上的点
ddata->r, ddata->c, curv_thr ) )
{
cvSeqPush( features, feat ); //是特征点进入特征点序列
}
else
free( ddata );
free( feat );
}
}
return features;
}
4.1
寻找极值点
每一个采样点都要与它相邻的像素点比较，看是否比它在图像域或尺度域的所有点的值大或者小。与它同尺度的相邻像素点有8个，上下相邻尺度的点共有2×9=18，共有26个像素点。也就在一个3×3的立方体内进行。搜索的过程是第二层开始到倒数第二层结束，共检测了octave组，每组S层。
/*
Determines whether a pixel is a scale-space extremum by comparing it to it's
3x3x3 pixel neighborhood.
*/
static int is_extremum( IplImage*** dog_pyr, int octv, int intvl, int r, int c )
{
double val = pixval32f( dog_pyr[octv][intvl], r, c );
int i, j, k;
/* check for maximum */
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;
}
/* check for minimum */
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;
}
4.2
准确定位特征点
以上的极值点搜索是在离散空间进行的，极值点不真正意义上的极值点。通过对空间尺度函数拟合，可以得到亚像素级像素点坐标。
尺度空间的Taylor展开式：
，其中
求导并令其为0，得到亚像素级：
对应的函数值为：
是一个三维矢量，矢量在任何一个方向上的偏移量大于0.5时，意味着已经偏离了原像素点，这样的特征坐标位置需要更新或者继续插值计算。算法实现过程中，为了保证插值能够收敛，设置了最大插值次数（lowe 设置了5次）。同时当
时（本文阈值采用了0.04/S）
，特征点才被保留，因为响应值过小的点，容易受噪声的干扰而不稳定。
对离散空间进行函数拟合（插值）：
/*
Performs one step of extremum interpolation. Based on Eqn. (3) in Lowe's
paper.
r,c 为特征点位置，xi,xr,xc,保存三个方向的偏移量
*/
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 };
dD = deriv_3D( dog_pyr, octv, intvl, r, c ); //计算三个方向的梯度
H = hessian_3D( dog_pyr, octv, intvl, r, c ); // 计算3维空间的hessian矩阵
H_inv = cvCreateMat( 3, 3, CV_64FC1 );
cvInvert( H, H_inv, CV_SVD ); //计算逆矩阵
cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );
cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 ); //广义乘法
cvReleaseMat( &dD );
cvReleaseMat( &H );
cvReleaseMat( &H_inv );
*xi = x[2];
*xr = x[1];
*xc = x[0];
}
/*
Interpolates a scale-space extremum's location and scale to subpixel
accuracy to form an image feature.
*/
static struct feature* interp_extremum( IplImage*** dog_pyr, int octv, //通过拟合求取准确的特征点位置
int intvl, int r, int c, int intvls,
double contr_thr )
{
struct feature* feat;
struct detection_data* ddata;
double xi, xr, xc, contr;
int i = 0;
while( i < SIFT_MAX_INTERP_STEPS ) //在最大迭代次数范围内进行
{
interp_step( dog_pyr, octv, intvl, r, c, &xi, &xr, &xc ); //插值后得到的三个方向的偏移量(xi,xr,xc)
if( ABS( xi ) < 0.5 && ABS( xr ) < 0.5 && ABS( xc ) < 0.5 )
break;
c += cvRound( xc ); //更新位置
r += cvRound( xr );
intvl += cvRound( xi );
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++;
}
/* ensure convergence of interpolation */
if( i >= SIFT_MAX_INTERP_STEPS )
return NULL;
contr = interp_contr( dog_pyr, octv, intvl, r, c, xi, xr, xc ); //计算插值后对应的函数值
if( ABS( contr ) < contr_thr / intvls ) //小于阈值（0.04/S）的点，则丢弃
return NULL;
feat = new_feature();
ddata = feat_detection_data( feat );
feat->img_pt.x = feat->x = ( c + xc ) * pow( 2.0, 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; // 层方向的偏移量
return feat;
}
待续...
5.计算不同尺度下特征点对应原图的位置
6.为特征点分配方向角
7.构建特征描述子
版权声明：本文为博主原创文章，未经博主允许不得转载。
时间： 2024-08-01 04:23:05