版權說明：本文檔由用戶提供并上傳，收益歸屬內(nèi)容提供方，若內(nèi)容存在侵權，請進行舉報或認領
文檔簡介
RobHess的SIFT算法實現(xiàn)理解及注釋
SIFT算法不用我多解釋了，這是一個很強大的算法，主要用于圖像配準和物體識別等領域，但是其計算量相比也比較大，性價比比較高的算法包括PCA-SIFT和SURF其中OpenCV提供了SURF算法，但是為了方便理解。這里給出了RobHess所實現(xiàn)的SIFT算法的實現(xiàn)以及注釋，結(jié)合我自己的理解，如果，您有關于SIFT算法不理解的地方咱們可以一起交流一下。或者您認為不詳細的地方提出來。
SIFT算法的主要實現(xiàn)在sift.c這個文件，其主要流程為：(1)首先創(chuàng)建初始圖像，即通過將圖像轉(zhuǎn)換為32位的灰度圖，然后將圖像使用三次插值來方大，之后通過高斯模糊處理(2)在此基礎上進行高斯金字塔的構(gòu)建以及高斯差分金字塔的構(gòu)建(3)對圖像進行極值點檢測(4)計算特征向量的尺度(5)調(diào)整圖像大小(6)計算特征的方向(7)計算描述子，其中包括計算二維方向直方圖并轉(zhuǎn)換直方圖為特征描述子首先給出sift算法的整體框架代碼：輸入?yún)?shù)：img為輸入圖像；feat為所要提取的特征指針；intvl指的是高斯金字塔和差分金字塔的層數(shù)；sigma指的是圖像初始化過程中高斯模糊所使用的參數(shù)；
contr_thr是歸一化之后的去除不穩(wěn)定特征的閾值；curv_thr指的是去除邊緣的特征的主曲率閾值；img_dbl是是否將圖像放大為之前的兩倍；descr_with用來計算特征描述子的方向直方圖的寬度；descr_hist_bins是直方圖中的條數(shù)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
*/
/*
構(gòu)建高斯尺度空間金字塔，頂層最小的為4像素
*/
init_img
=
create_init_img(
img,
img_dbl,
sigma
);
octvs
=
log(
double
MIN(
init_img->width,
init_img->height
)
)
/
log(2.0)
-
2;
//構(gòu)建高斯金字塔和高斯差分金字塔
gauss_pyr
=
build_gauss_pyr(
init_img,
octvs,
intvls,
sigma
);
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
);
//畫出去除低對比度的極值點
//draw_extrempoint(img
,
features);
//計算特征向量的尺度
calc_feature_scales(
features,
sigma,
intvls
);
if(
img_dbl
)
adjust_for_img_dbl(
features
);
//計算特征的方向
calc_feature_oris(
features,
gauss_pyr
);
//計算描述子，包括計算二維方向直方圖和轉(zhuǎn)換其為特征描述子
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
=
static_cast<feature
*>(
calloc(
n,
sizeof(struct
feature)
)
);
*feat
=
static_cast<feature
*>(
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)初始化圖像輸入?yún)?shù)：這里不需要解釋了該函數(shù)主要用來初始化圖像，轉(zhuǎn)換圖像為32位灰度圖以及進行高斯模糊。static
IplImage*
create_init_img(
IplImage*
img,
int
img_dbl,
double
sigma
)
{
IplImage*
gray,
*
dbl;
float
sig_diff;
gray
=
convert_to_gray32(
img
);
if(
img_dbl
)
{
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
);
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;
}
}
(2)構(gòu)建高斯金字塔輸入?yún)?shù)：octvs是高斯金字塔的組invls是高斯金字塔的層數(shù)sigma是初始的高斯模糊參數(shù)，后續(xù)也通過它計算每一層所使用的sigma<span
style="font-size:13px;">static
IplImage***
build_gauss_pyr(
IplImage*
base,
int
octvs,int
intvls,
double
sigma
)
{
IplImage***
gauss_pyr;
double*
sig
=
static_cast<double
*>(
calloc(
intvls
+
3,
sizeof(double))
);
double
sig_total,
sig_prev,
k;
int
i,
o;
gauss_pyr
=
static_cast<IplImage
***>(
calloc(
octvs,
sizeof(
IplImage**
)
)
);
for(
i
=
0;
i
<
octvs;
i++
)
gauss_pyr[i]
=
static_cast<IplImage
**>(
calloc(
intvls
+
3,
sizeof(
IplImage*
)
)
);
/*
precompute
Gaussian
sigmas
using
the
following
formula:
預計算每次高斯模糊的sigma
\sigma_{total}^2
=
\sigma_{i}^2
+
\sigma_{i-1}^2
*/
sig[0]
=
sigma;
k
=
pow(
2.0,
1.0
/
intvls
);
for(
i
=
1;
i
<
intvls
+
3;
i++
)
{
sig_prev
=
pow(
k,
i
-
1
)
*
sigma;
sig_total
=
sig_prev
*
k;
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
)
gauss_pyr[o][i]
=
cvCloneImage(base);
/*
base
of
new
octvave
is
halved
image
from
end
of
previous
octave
*/
//每一組的第一層都是通過對前面一組的最上面一層的降采樣實現(xiàn)的
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
*/
//每一組的其他層則使通過使用不同sigma的高斯模糊來進行處理
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]
);
}
}
free(
sig
);
return
gauss_pyr;
}</span>
降采樣處理輸入?yún)?shù)：不解釋這就是降采樣，其實就是將圖像通過最近鄰算法縮小為原來的一半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;
}
(3)構(gòu)建高斯差分金字塔輸入?yún)?shù)：不解釋了參見上面的說明即可實際上差分金字塔的構(gòu)成是通過對相鄰層的圖像進行相減獲得的<span
style="font-size:16px;">static
IplImage***
build_dog_pyr(
IplImage***
gauss_pyr,
int
octvs,
int
intvls
)
{
IplImage***
dog_pyr;
int
i,
o;
dog_pyr
=
static_cast<IplImage
***>(
calloc(
octvs,
sizeof(
IplImage**
)
)
);
for(
i
=
0;
i
<
octvs;
i++
)
dog_pyr[i]
=
static_cast<IplImage
**>(
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
);
}
return
dog_pyr;
}</span>
(4)極值點檢測輸入?yún)?shù)：contr_thr是去除對比度低的點所采用的閾值curv_thr是去除邊緣特征的閾值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++
)
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
)
if(
is_extremum(
dog_pyr,
o,
i,
r,
c
)
)
{
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],
ddata->r,
ddata->c,
curv_thr
)
)
{
cvSeqPush(
features,
feat
);
}
else
free(
ddata
);
free(
feat
);
}
}
return
features;
}
SIFT_IMG_BORDER是預定義的圖像邊緣；通過和對比度閾值比較去掉低對比度的點；而通過is_extremum來判斷是否為極值點，如果是則通過極值點插值的方式獲取亞像素的極值點的位置。然后通過is_too_eage_like和所給的主曲率閾值判斷是否為邊緣點*判斷是否為極值點其原理為:通過和高斯金字塔的上一層的9個像素+本層的除了本像素自己的其他的8個像素和下一層的9個像素進行比較看是否為這26個像素中最小的一個或者是否為最大的一個，如果是則為極值點。static
int
is_extremum(
IplImage***
dog_pyr,
int
octv,
int
intvl,
int
r,
int
c
)
{
float
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;
}
*獲取亞像素的極值點的位置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;//分別為亞像素的intval,row,col的偏移offset，和對比度
int
i
=
0;
while(
i
<
SIFT_MAX_INTERP_STEPS
)//重新確定極值點并重新定位的操作只能循環(huán)
5次
{
interp_step(
dog_pyr,
octv,
intvl,
r,
c,
&xi,
&xr,
&xc
);
if(
ABS(
xi
)
<
0.5
&&
ABS(
xr
)
<
0.5
&&
ABS(
xc
)
<
0.5
)//如果滿足條件就停止尋找
break;
//否則繼續(xù)尋找極值點
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++;
}
//確保極值點是經(jīng)過最大5步找到的
/*
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
)
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;
}
*獲取亞像素位置中所用到的函數(shù)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
};
//計算三維偏導數(shù)
dD
=
deriv_3D(
dog_pyr,
octv,
intvl,
r,
c
);
//計算三維海森矩陣
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
);
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];
}
*計算三維偏導數(shù)計算在x和y方向上的偏導數(shù)，高斯差分尺度空間金字塔中像素的尺度實際上在離散數(shù)據(jù)中計算偏導數(shù)是通過相鄰像素的相減來計算的比如說計算x方向的偏導數(shù)dx，則通過該向所的x方向的后一個減去前一個然后除以2即可求的dxstatic
CvMat*
deriv_3D(
IplImage***
dog_pyr,
int
octv,
int
intvl,
int
r,
int
c
)
{
CvMat*
dI;
double
dx,
dy,
ds;
dx
=
(
pixval32f(
dog_pyr[octv][intvl],
r,
c+1
)
-
pixval32f(
dog_pyr[octv][intvl],
r,
c-1
)
)
/
2.0;
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;
}
*計算三維海森矩陣不需要講什么，其實就是計算二次導數(shù)，計算方法也和一次導數(shù)的計算如出一轍。然后將結(jié)果放入到一個矩陣中去。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
=
(
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;
}
*計算插入像素的對比度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
};
cvInitMatHeader(
&X,
3,
1,
CV_64FC1,
x,
CV_AUTOSTEP
);
cvInitMatHeader(
&T,
1,
1,
CV_64FC1,
t,
CV_AUTOSTEP
);
dD
=
deriv_3D(
dog_pyr,
octv,
intvl,
r,
c
);
cvGEMM(
dD,
&X,
1,
NULL,
0,
&T,
CV_GEMM_A_T
);
cvReleaseMat(
&dD
);
return
pixval32f(
dog_pyr[octv][intvl],
r,
c
)
+
t[0]
*
0.5;
}
其中cvGEMM是矩陣的通用計算函數(shù)，至于CV_GEMM_A_T是計算dD的轉(zhuǎn)置矩陣放入T中*去除邊緣相應
通過計算所在特征向量的主曲率半徑來判斷特征是邊緣的從而導致不穩(wěn)定
即去除邊緣響應static
int
is_too_edge_like(
IplImage*
dog_img,
int
r,
int
c,
int
curv_thr
)
{
double
d,
dxx,
dyy,
dxy,
tr,
det;
/*
principal
curvatures
are
computed
using
the
trace
and
det
of
Hessian
*/
d
=
pixval32f(dog_img,
r,
c);
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;
if(
tr
*
tr
/
det
<
(
curv_thr
+
1.0
)*(
curv_thr
+
1.0
)
/
curv_thr
)
return
0;
return
1;
}
(4)計算特征向量的尺度實際上是通過最初的sigma來獲得每一層每一組的尺度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++
)
{
feat
=
CV_GET_SEQ_ELEM(
struct
feature,
features,
i
);
ddata
=
feat_detection_data(
feat
);
intvl
=
ddata->intvl
+
ddata->subintvl;
feat->scl
=
sigma
*
pow(
2.0,
ddata->octv
+
intvl
/
intvls
);
ddata->scl_octv
=
sigma
*
pow(
2.0,
intvl
/
intvls
);
}
}
(5)調(diào)整圖像特征坐標、尺度、點的坐標大小為原來的一半static
void
adjust_for_img_dbl(
CvSeq*
features
)
{
struct
feature*
feat;
int
i,
n;
n
=
features->total;
for(
i
=
0;
i
<
n;
i++
)
{
feat
=
CV_GET_SEQ_ELEM(
struct
feature,
features,
i
);
feat->x
/=
2.0;
feat->y
/=
2.0;
feat->scl
/=
2.0;
feat->img_pt.x
/=
2.0;
feat->img_pt.y
/=
2.0;
}
}
(6)給每一個圖像特征向量計算規(guī)范化的方向static
void
calc_feature_oris(
CvSeq*
features,
IplImage***
gauss_pyr
)
{
struct
feature*
feat;
struct
detection_data*
ddata;
double*
hist;
double
omax;
int
i,
j,
n
=
features->total;
//遍歷整個檢測出來的特征點，計算每個特征點的直方圖，然后平滑直方圖去除突變，然后找到每一個特征點的主方向，并加入到好的方向特征數(shù)組中去
for(
i
=
0;
i
<
n;
i++
)
{
feat
=
static_cast<feature
*>(
malloc(
sizeof(
struct
feature
)
)
);
cvSeqPopFront(
features,
feat
);
ddata
=
feat_detection_data(
feat
);
//計算給定的某個像素的灰度方向直方圖
hist
=
ori_hist(
gauss_pyr[ddata->octv][ddata->intvl],
ddata->r,
ddata->c,
SIFT_ORI_HIST_BINS,
cvRound(
SIFT_ORI_RADIUS
*
ddata->scl_octv
),
SIFT_ORI_SIG_FCTR
*
ddata->scl_octv
);
for(
j
=
0;
j
<
SIFT_ORI_SMOOTH_PASSES;
j++
)
smooth_ori_hist(
hist,
SIFT_ORI_HIST_BINS
);
omax
=
dominant_ori(
hist,
SIFT_ORI_HIST_BINS
);
//描述子向量元素門限化
add_good_ori_features(
features,
hist,
SIFT_ORI_HIST_BINS,
omax
*
SIFT_ORI_PEAK_RATIO,
feat
);
free(
ddata
);
free(
feat
);
free(
hist
);
}
}
*對所給像素計算灰度方向直方圖
以關鍵點為中心的鄰域窗口內(nèi)采樣，并用直方圖統(tǒng)計鄰域像素的梯度
方向。梯度直方圖的范圍是0～360度，其中每10度一個柱，總共36個柱static
double*
ori_hist(
IplImage*
img,
int
r,
int
c,
int
n,
int
rad,
double
sigma)
{
double*
hist;
double
mag,
ori,
w,
exp_denom,
PI2
=
CV_PI
*
2.0;
int
bin,
i,
j;
hist
=
static_cast<double
*>(
calloc(
n,
sizeof(
double
)
)
);
exp_denom
=
2.0
*
sigma
*
sigma;
for(
i
=
-rad;
i
<=
rad;
i++
)
for(
j
=
-rad;
j
<=
rad;
j++
)
if(
calc_grad_mag_ori(
img,
r
+
i,
c
+
j,
&mag,
&ori
)
)
{
w
=
exp(
-(
i*i
+
j*j
)
/
exp_denom
);
bin
=
cvRound(
n
*
(
ori
+
CV_PI
)
/
PI2
);
bin
=
(
bin
<
n
)?
bin
:
0;
hist[bin]
+=
w
*
mag;
}
return
hist;
}
*計算所給像素的梯度大小和方向
每一個小格都代表了特征點鄰域所在的尺度空間的一個像素，箭頭方向代表了像素梯
度方向，箭頭長度代表該像素的幅值也就是梯度的值static
int
calc_grad_mag_ori(
IplImage*
img,
int
r,
int
c,
double*
mag,
double*
ori
)
{
double
dx,
dy;
if(
r
>
0
&&
r
<
img->height
-
1
&&
c
>
0
&&
c
<
img->width
-
1
)
{
dx
=
pixval32f(
img,
r,
c+1
)
-
pixval32f(
img,
r,
c-1
);
dy
=
pixval32f(
img,
r-1,
c
)
-
pixval32f(
img,
r+1,
c
);
*mag
=
sqrt(
dx*dx
+
dy*dy
);
*ori
=
atan2(
dy,
dx
);
return
1;
}
else
return
0;
}
*對方向直方圖進行高斯模糊
使用高斯函數(shù)對直方圖進行平滑，減少突變的影響。static
void
smooth_ori_hist(
double*
hist,
int
n
)
{
double
prev,
tmp,
h0
=
hist[0];
int
i;
prev
=
hist[n-1];
for(
i
=
0;
i
<
n;
i++
)
{
tmp
=
hist[i];
hist[i]
=
0.25
*
prev
+
0.5
*
hist[i]
+
0.25
*
(
(
i+1
==
n
)?
h0
:
hist[i+1]
);
prev
=
tmp;
}
}
*在直方圖中找到主方向的梯度
利用關鍵點鄰域像素的梯度方向分布特性為每個關鍵點指定方向參數(shù)，使算子具備
旋轉(zhuǎn)不變性。static
double
dominant_ori(
double*
hist,
int
n
)
{
double
omax;
int
maxbin,
i;
omax
=
hist[0];
maxbin
=
0;
for(
i
=
1;
i
<
n;
i++
)
if(
hist[i]
>
omax
)
{
omax
=
hist[i];
maxbin
=
i;
}
return
omax;
}
*將大于某一個梯度大小閾值的特征向量加入到直方圖中去n為方向的個數(shù)<span
style="font-size:18px;">mag_thr描述子向量門限一般取0.2</span>
static
void
add_good_ori_features(
CvSeq*
features,
double*
hist,
int
n,
double
mag_thr,
struct
feature*
feat
)
{
struct
feature*
new_feat;
double
bin,
PI2
=
CV_PI
*
2.0;
int
l,
r,
i;
for(
i
=
0;
i
<
n;
i++
)
{
l
=
(
i
==
0
)?
n
-
1
:
i-1;
r
=
(
i
+
1
)
%
n;
//描述子向量門限化，一般門限取0.2
if(
hist[i]
>
hist[l]
&&
hist[i]
>
hist[r]
&&
hist[i]
>=
mag_thr
)
{
bin
=
i
+
interp_hist_peak(
hist[l],
hist[i],
hist[r]
);
bin
=
(
bin
<
0
)?
n
+
bin
:
(
bin
>=
n
)?
bin
-
n
:
bin;
new_feat
=
clone_feature(
feat
);
new_feat->ori
=
(
(
PI2
*
bin
)
/
n
)
-
CV_PI;
cvSeqPush(
features,
new_feat
);
free(
new_feat
);
}
}
}
(7)計算特征描述子static
void
compute_descriptors(
CvSeq*
features,
IplImage***
gauss_pyr,
int
d,
int
n)
{
struct
feature*
feat;
struct
detection_data*
ddata;
double***
hist;
int
i,
k
=
features->total;
for(
i
=
0;
i
<
k;
i++
)
{
feat
=
CV_GET_SEQ_ELEM(
struct
feature,
features,
i
);
ddata
=
feat_detection_data(
feat
);
//計算二維方向直方圖
hist
=
descr_hist(
gauss_pyr[ddata->octv][ddata->intvl],
ddata->r,
ddata->c,
feat->ori,
ddata->scl_octv,
d,
n
);
//將二維方向直方圖轉(zhuǎn)換為特征描述子
hist_to_descr(
hist,
d,
n,
feat
);
release_descr_hist(
&hist,
d
);
}
}
*計算二維方向直方圖static
double***
descr_hist(
IplImage*
img,
int
r,
int
c,
double
ori,
double
scl,
int
d,
int
n
)
{
double***
hist;
double
cos_t,
sin_t,
hist_width,
exp_denom,
r_rot,
c_rot,
grad_mag,
grad_ori,
w,
rbin,
cbin,
obin,
bins_per_rad,
PI2
=
2.0
*
CV_PI;
int
radius,
i,
j;
hist
=
static_cast<double
***>(
calloc(
d,
sizeof(
double**
)
)
);
for(
i
=
0;
i
<
d;
i++
)
{
hist[i]
=static_cast<double
**>(
calloc(
d,
sizeof(
double*
)
)
);
for(
j
=
0;
j
<
d;
j++
)
hist[i][j]
=
static_cast<double
*>(
calloc(
n,
sizeof(
double
)
)
);
}
cos_t
=
cos(
ori
);
sin_t
=
sin(
ori
);
bins_per_rad
=
n
/
PI2;
exp_denom
=
d
*
d
*
0.5;
hist_width
=
SIFT_DESCR_SCL_FCTR
*
scl;
radius
=
hist_width
*
sqrt(2.0)
*
(
d
+
1.0
)
*
0.5
+
0.5;
for(
i
=
-radius;
i
<=
radius;
i++
)
for(
j
=
-radius;
j
<=
radius;
j++
)
{
/*
即將坐標移至關鍵點主方向
計算采用的直方圖數(shù)組中相對于方向旋轉(zhuǎn)的坐標
Calculate
sample's
histogram
array
coords
rotated
relative
to
ori.
Subtract
0.5
so
samples
that
fall
e.g.
in
the
center
of
row
1
(i.e.
r_rot
=
1.5)
have
full
weight
placed
in
row
1
after
interpolation.
*/
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
)
if(
calc_grad_mag_ori(
img,
r
+
i,
c
+
j,
&grad_mag,
&grad_ori
))
{
grad_ori
-=
ori;
while(
grad_ori
<
0.0
)
grad_ori
+=
PI2;
while(
grad_ori
>=
PI2
)
grad_ori
-=
PI2;
obin
=
grad_ori
*
bins_per_rad;
w
=
exp(
-(c_rot
*
c_rot
+
r_rot
*
r_rot)
/
exp_denom
);
interp_hist_entry(
hist,
rbin,
cbin,
obin,
grad_mag
*
w,
d,
n
);
}
}
return
hist;
}
*插入一個entry進入到方向直方圖中從而形成特征描述子這個，我也不怎么明白。。。static
void
interp_hist_entry(
double***
hist,
double
rbin,
double
cbin,
double
obin,
double
mag,
int
d,
int
n
)
{
double
d_r,
d_c,
d_o,
v_r,
v_c,
v_o;
double**
row,
*
h;
int
r0,
c0,
o0,
rb,
cb,
ob,
r,
c,
o;
r0
=
cvFloor(
rbin
);
c0
=
cvFloor(
cbin
);
o0
=
cvFloor(
obin
);
d_r
=
rbin
-
r0;
d_c
=
cbin
-
c0;
d_o
=
obin
-
o0;
/*
The
entry
is
distributed
into
up
to
8
bins.
Each
entry
into
a
bin
is
multiplied
by
a
weight
of
1
-
d
for
each
dimension,
where
d
is
the
distance
from
the
center
value
of
the
bin
measured
in
bin
units.
*/
for(
r
=
0;
r
<=
1;
r++
)
{
rb
=
r0
+
r;
if(
rb
>=
0
&&
rb
<
d
)
{
v_r
=
mag
*
(
(
r
==
0
)?
1.0
-
d_r
:
d_r
);
row
=
hist[rb];
for(
c
=
0;
c
<=
1;
c++
)
{
cb
=
c0
+
c;
if(
cb
>=
0
&&
cb
<
d
)
{
v_c
=
v_r
*
(
(
c
==
0
)?
1.0
-
d_c
:
d_c
);
h
=
row[cb];
for(
o
=
0;
o
<=
1;
o++
)
{
ob
=
(
o0
+
o
)
%
n;
v_o
=
v_c
*
(
(
o
==
0
)?
1.0
-
d_o
:
d_o
);
溫馨提示
- 1. 本站所有資源如無特殊說明，都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
- 2. 本站的文檔不包含任何第三方提供的附件圖紙等，如果需要附件，請聯(lián)系上傳者。文件的所有權益歸上傳用戶所有。
- 3. 本站RAR壓縮包中若帶圖紙，網(wǎng)頁內(nèi)容里面會有圖紙預覽，若沒有圖紙預覽就沒有圖紙。
- 4. 未經(jīng)權益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
- 5. 人人文庫網(wǎng)僅提供信息存儲空間，僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護處理，對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯，并不能對任何下載內(nèi)容負責。
- 6. 下載文件中如有侵權或不適當內(nèi)容，請與我們聯(lián)系，我們立即糾正。
- 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。
最新文檔
- 高效溝通協(xié)作機制建立方案
- 鄉(xiāng)村環(huán)境綜合整治技術作業(yè)指導書
- 電力行業(yè)供電安全告知書
- 房屋買賣按揭合同
- 商業(yè)場所租賃使用協(xié)議及設備設施管理細則協(xié)議
- 智能辦公系統(tǒng)集成方案簽署協(xié)議
- 高考語文復習-文言文重點字詞解析練習
- 高考英語整句翻譯漢譯英專題訓練500題（含答案）
- 新品手機使用說明手冊
- 企業(yè)研發(fā)創(chuàng)新基金合作協(xié)議
- 2024年鄭州信息科技職業(yè)學院高職單招（英語／數(shù)學／語文）筆試歷年參考題庫含答案解析
- 藍牙基礎知識全解課件
- 運動損傷預防與處理的案例分析
- 第四次工業(yè)革命課件
- 2023-2024學年西安市高二數(shù)學第一學期期末考試卷附答案解析
- 企業(yè)2024年年度安全教育培訓計劃
- 《微生物限度檢查法》課件
- Project-培訓教學課件
- 秋風詞賞析課件古詩詞賞析
- 福特F-150猛禽說明書
- DB3402-T 59-2023 露天礦山無人駕駛礦車作業(yè)通用要求
評論
0/150
提交評論