sift演算法c語言實現
阿新 • • 發佈：2019-02-17
前段時間在做三維測量方面的研究，需要得到物體表面三維資料，sift演算法是立體匹配中的經典演算法，下面是對RobHess的SIFT原始碼的註釋。部分內容參考網上，在這裡向各位大神表示感謝！
http://blog.csdn.net/lsh_2013/article/details/46826141
標頭檔案及函式宣告
#include "sift.h"
#include "imgfeatures.h"
#include "utils.h"
#include <cxcore.h>
#include <cv.h>
//將原圖轉換為32位灰度圖並歸一化，然後進行一次高斯平滑，並根據引數img_dbl決定 函式實現
主函式
/*使用使用者指定的引數在影象中提取SIFT特徵點
引數：
img：輸入影象
feat：儲存特徵點的陣列的指標，此陣列的記憶體將在本函式中被分配，使用完後必須在調用出釋放：free(*feat)
intvls：每組的層數
sigma：初始高斯平滑引數σ
contr_thr：對比度閾值，針對歸一化後的影象，用來去除不穩定特徵
curv_thr：去除邊緣的特徵的主曲率閾值
img_dbl：是否將影象放大為之前的兩倍
descr_width：特徵描述過程中，計算方向直方圖時，將特徵點附近劃分為descr_width*descr_width個區域，每個區域生成一個直方圖
descr_hist_bins：特徵描述過程中，每個直方圖中bin的個數
返回值：提取的特徵點個數，若返回-1表明提取失敗*/
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;//三級指標，高斯金字塔影象組，DoG金字塔影象組
CvMemStorage* storage;//儲存器
CvSeq* features;//儲存特徵點的序列，序列中存放的是struct feature型別的指標
int octvs, i, n = 0;
//輸入引數檢查
if( ! img )
fatal_error( "NULL pointer error, %s, line %d", __FILE__, __LINE__ );
if( ! feat )
fatal_error( "NULL pointer error, %s, line %d", __FILE__, __LINE__ );
/*步驟一：建立尺度空間，即建立高斯差分(DoG)金字塔dog_pyr
將原圖轉換為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 );
//建立高斯差分(DoG)金字塔dog_pyr，是一個octvs*(intvls+2)的影象陣列
dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );
/*步驟二：在尺度空間中檢測極值點，並進行精確定位和篩選建立預設大小的記憶體儲存器*/
storage = cvCreateMemStorage( 0 );
//在尺度空間中檢測極值點，通過插值精確定位，去除低對比度的點，去除邊緣點，
//返回檢測到的特徵點序列
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倍時，特徵點檢測完之後呼叫)
/*步驟三：特徵點方向賦值，完成此步驟後，每個特徵點有三個資訊：位置、尺度、方向*/
//計算每個特徵點的梯度直方圖，找出其主方向，若一個特徵點有不止一個主方向，將其分為
//兩個特徵點
calc_feature_oris( features, gauss_pyr );
/*步驟四：計算特徵描述子*/
//計算特徵點序列中每個特徵點的特徵描述子向量
compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins );
//按特徵點尺度的降序排列序列中的元素的順序，feature_cmp是自定義的比較函式
cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL );
//將CvSeq型別的特徵點序列features轉換為通用的struct feature型別的陣列feat
n = features->total;//特徵點個數
*feat = calloc( n, sizeof(struct feature) );//分配控制元件
//將序列features中的元素拷貝到陣列feat中，返回陣列指標給feat
*feat = cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ );
//釋放特徵點陣列feat中所有特徵點的feature_data成員，因為此成員中的資料在檢測完特徵
//點後就沒用了
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;
}尺度空間構造
/*將原圖轉換為32位灰度圖並歸一化，然後進行一次高斯平滑，並根據引數img_dbl決定是否將
影象尺寸放大為原圖的2倍
引數：
img：輸入的原影象
img_dbl：是否將影象放大為之前的兩倍
sigma：初始高斯平滑引數σ
返回值：初始化完成的影象*/
static IplImage* create_init_img( IplImage* img, int img_dbl, double sigma )
{
IplImage* gray, * dbl;
float sig_diff;
//呼叫函式，將輸入影象轉換為32位灰度圖，並歸一化
gray = convert_to_gray32( img );
if( img_dbl ) //若設定了將影象放大為原圖的2倍
{
//將影象長寬擴充套件一倍時，便有了底-1層，該層尺度為：
sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4 );
//建立放大影象
dbl = cvCreateImage( cvSize( img->width*2, img->height*2 ),
IPL_DEPTH_32F, 1 );
//放大原圖的尺寸
cvResize( gray, dbl, CV_INTER_CUBIC );
//高斯平滑，高斯核在x,y方向上的標準差都是sig_diff
cvSmooth( dbl, dbl, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );
cvReleaseImage( &gray );
return dbl;
}
else//不用放大為原圖的2倍
{
//計算第0層的尺度
sig_diff = sqrt( sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA );
//高斯平滑，高斯核在x,y方向上的標準差都是sig_diff
cvSmooth( gray, gray, CV_GAUSSIAN, 0, 0, sig_diff, sig_diff );
return gray;
}
}
/*將輸入影象轉換為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 )//若原圖本身就是單通道，直接克隆原圖
gray8 = cvClone( img );
else//若原圖是3通道影象
{
gray8 = cvCreateImage( cvGetSize(img), IPL_DEPTH_8U, 1 );//建立8位單通道影象
cvCvtColor( img, gray8, CV_BGR2GRAY );//將原圖轉換為8為單通道影象
}
//然後將8為單通道影象gray8轉換為32位單通道影象，並進行歸一化處理(除以255)
cvConvertScale( gray8, gray32, 1.0 / 255.0, 0 );
//釋放臨時影象
cvReleaseImage( &gray8 );
//返回32位單通道影象
return gray32;
}
/*根據輸入引數建立高斯金字塔
引數：
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個元素，每個元素是一組影象的首指標
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;//返回高斯金字塔
}
/*對輸入影象做下采樣生成其四分之一大小的影象(每個維度上減半)，使用最近鄰差值方法
引數：
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 );//尺寸變換
return smaller;
}
/*通過對高斯金字塔中每相鄰兩層影象相減來建立高斯差分金字塔
引數：
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 );
}
return dog_pyr;//返回高斯差分金字塔
}
區域性空間極值點檢測
/*在尺度空間中檢測極值點，通過插值精確定位，去除低對比度的點，去除邊緣點，返回檢測到
的特徵點序列
引數：
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 )
{
CvSeq* features;//特徵點序列
double prelim_contr_thr = 0.5 * contr_thr / intvls;//畫素的對比度閾值
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;//返回特徵點序列
}
/*通過在尺度空間中將一個畫素點的值與其周圍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;
}剔除不穩定點，精確定位關鍵點位置
/*通過亞畫素級插值進行極值點精確定位(修正極值點座標)，並去除低對比度的極值點，
將修正後的特徵點組成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 )
{
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 );
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;//返回特徵點指標
}
/*進行一次極值點差值，計算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方向上偏移量
}
/*在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;
}
/*在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;
}
/*計算被插值點的對比度：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;
}
/*為一個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;
}
/*去除邊緣響應，即通過計算主曲率比值判斷某點是否邊緣點
引數：
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;//是邊緣點
}
/*計算特徵點序列中每個特徵點的尺度
引數：
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 );