版权说明：本文档由用户提供并上传，收益归属内容提供方，若内容存在侵权，请进行举报或认领
文档简介
RobHess的SIFT算法实现理解及注释
SIFT算法不用我多解释了，这是一个很强大的算法，主要用于图像配准和物体识别等领域，但是其计算量相比也比较大，性价比比较高的算法包括PCA-SIFT和SURF其中OpenCV提供了SURF算法，但是为了方便理解。这里给出了RobHess所实现的SIFT算法的实现以及注释，结合我自己的理解，如果，您有关于SIFT算法不理解的地方咱们可以一起交流一下。或者您认为不详细的地方提出来。
SIFT算法的主要实现在sift.c这个文件，其主要流程为：(1)首先创建初始图像，即通过将图像转换为32位的灰度图，然后将图像使用三次插值来方大，之后通过高斯模糊处理(2)在此基础上进行高斯金字塔的构建以及高斯差分金字塔的构建(3)对图像进行极值点检测(4)计算特征向量的尺度(5)调整图像大小(6)计算特征的方向(7)计算描述子，其中包括计算二维方向直方图并转换直方图为特征描述子首先给出sift算法的整体框架代码：输入参数：img为输入图像；feat为所要提取的特征指针；intvl指的是高斯金字塔和差分金字塔的层数；sigma指的是图像初始化过程中高斯模糊所使用的参数；
contr_thr是归一化之后的去除不稳定特征的阈值；curv_thr指的是去除边缘的特征的主曲率阈值；img_dbl是是否将图像放大为之前的两倍；descr_with用来计算特征描述子的方向直方图的宽度；descr_hist_bins是直方图中的条数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
*/
/*
构建高斯尺度空间金字塔，顶层最小的为4像素
*/
init_img
=
create_init_img(
img,
img_dbl,
sigma
);
octvs
=
log(
double
MIN(
init_img->width,
init_img->height
)
)
/
log(2.0)
-
2;
//构建高斯金字塔和高斯差分金字塔
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
);
//计算描述子，包括计算二维方向直方图和转换其为特征描述子
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)初始化图像输入参数：这里不需要解释了该函数主要用来初始化图像，转换图像为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)构建高斯金字塔输入参数：octvs是高斯金字塔的组invls是高斯金字塔的层数sigma是初始的高斯模糊参数，后续也通过它计算每一层所使用的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
*/
//每一组的第一层都是通过对前面一组的最上面一层的降采样实现的
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>
降采样处理输入参数：不解释这就是降采样，其实就是将图像通过最近邻算法缩小为原来的一半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)构建高斯差分金字塔输入参数：不解释了参见上面的说明即可实际上差分金字塔的构成是通过对相邻层的图像进行相减获得的<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)极值点检测输入参数：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
)//重新确定极值点并重新定位的操作只能循环
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;
//否则继续寻找极值点
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++;
}
//确保极值点是经过最大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;
}
*获取亚像素位置中所用到的函数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
);
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];
}
*计算三维偏导数计算在x和y方向上的偏导数，高斯差分尺度空间金字塔中像素的尺度实际上在离散数据中计算偏导数是通过相邻像素的相减来计算的比如说计算x方向的偏导数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;
}
*计算三维海森矩阵不需要讲什么，其实就是计算二次导数，计算方法也和一次导数的计算如出一辙。然后将结果放入到一个矩阵中去。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是矩阵的通用计算函数，至于CV_GEMM_A_T是计算dD的转置矩阵放入T中*去除边缘相应
通过计算所在特征向量的主曲率半径来判断特征是边缘的从而导致不稳定
即去除边缘响应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)调整图像特征坐标、尺度、点的坐标大小为原来的一半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)给每一个图像特征向量计算规范化的方向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;
//遍历整个检测出来的特征点，计算每个特征点的直方图，然后平滑直方图去除突变，然后找到每一个特征点的主方向，并加入到好的方向特征数组中去
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
);
}
}
*对所给像素计算灰度方向直方图
以关键点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度
方向。梯度直方图的范围是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;
}
*对方向直方图进行高斯模糊
使用高斯函数对直方图进行平滑，减少突变的影响。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;
}
}
*在直方图中找到主方向的梯度
利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备
旋转不变性。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为方向的个数<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
);
//将二维方向直方图转换为特征描述子
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++
)
{
/*
即将坐标移至关键点主方向
计算采用的直方图数组中相对于方向旋转的坐标
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. 本站的文档不包含任何第三方提供的附件图纸等，如果需要附件，请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸，网页内容里面会有图纸预览，若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对用户上传分享的文档内容本身不做任何修改或编辑，并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容，请与我们联系，我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 凿井勘查合同范例
- 劳务损伤赔偿合同范本
- 化工生产合同范本
- 2024年中国动漫博物馆（杭州）招聘考试真题
- 2024年重庆永川区五间镇招聘公益性岗位人员笔试真题
- 乡下房屋转卖合同范本
- gf分包合同范本
- 修路合同范本简版
- 出售小区公共用地合同范本
- 北京三室一厅租房合同范本
- 安全管理工作中形式主义及防止对策
- 2024年郑州信息科技职业学院高职单招（英语／数学／语文）笔试历年参考题库含答案解析
- 2023-2024学年西安市高二数学第一学期期末考试卷附答案解析
- 学校保密教育培训课件
- 班组文化是企业文化建设的核心
- Project-培训教学课件
- 福建省服务区标准化设计指南
- 秋风词赏析课件古诗词赏析
- 销售人员薪酬设计实例 薪酬制度设计 薪酬设计方案 设计案例全套
- 福特F-150猛禽说明书
- 征地搬迁基本要求及工作技巧课件
评论
0/150
提交评论