VLFeat是一个很好用的开源库，其中实现了计算机视觉常用的算法，使用的语言是C和matlab。
官网：http://www.vlfeat.org/index.html
在官网下载最新版本后，在matlab中添加路径即可。
1，Dense Sift
在文章《sift特征提取算法》中提到，提取图像的sift特征分4步：构建DOG尺度空间；关键点定位；关键点方向赋值；生成描述子。
这里产生的sift特征点是sparse sift，而在实际实验中，使用较多的则是dense sift。
dense sift省去了前3步，即直接指定关键点位置和描述子采样区域，计算sift特征。
主要过程是：
1，用一个patch在图像上以一定步长step滑动，代码中step=1，这个patch就是描述子采样区域，patch size是4bins*4bins，bin size可以自己指定，代码中是3pixels*3pixels。这里说的bin对应到《sift特征提取》中的第4步就是指子区域area。图中的bounding box是sift特征点的范围。
2，计算每个像素点的梯度（同sparse sift），统计每个bin内的像素点在8个方向上的梯度直方图，这样就生成了4*4*8维的sift特征。
在matlab中直接调用vl_dsift:
% 读入图像
I = vl_impattern(‘roofs1‘) ;
I = single(vl_imdown(rgb2gray(I))) ;
binSize = 8 ;
magnif = 3 ;
% 得到指定尺度的高斯图像
Is = vl_imsmooth(I, sqrt((binSize/magnif)^2 - .25)) ;
% 计算dense sift
[frame, descr] = vl_dsift(Is, ‘size‘, binSize) ;
% frame中存储每个sift点的坐标，descr存储每个sift点的128维特征向量主要实现代码在vl_dsift.c, dsift.c中
参数说明：
size: 是每个bin的大小，代码中是3pixels*3pixels
step: patch移动步长
bounds: sift点的区域，含minX, minY, maxX,maxY四个值
norm: 本来是存在frame中的一列，存储descriptor的值的和，用于normalization
geometry: [4 4 8]，在4bins*4bins的范围内，对每个bin统计8个方向的梯度直方图
verbose: 一般bin大小和patch大小都取正方形，也就是长和宽一致，如果打开这个verbose开关，则要对sizeX和sizeY分别赋值。
自定义数据结构在文件dsift.h中：
typedef struct VlDsiftKeypoint_
{
double x ; /**< x coordinate */
double y ; /**< y coordinate */
double s ; /**< scale */
double norm ; /**< SIFT descriptor norm */
} VlDsiftKeypoint ;
/** @brief Dense SIFT descriptor geometry */
typedef struct VlDsiftDescriptorGeometry_
{
int numBinT ; /**< number of orientation bins */
int numBinX ; /**< number of bins along X */
int numBinY ; /**< number of bins along Y */
int binSizeX ; /**< size of bins along X */
int binSizeY ; /**< size of bins along Y */
} VlDsiftDescriptorGeometry ;
typedef struct VlDsiftFilter_
{
int imWidth ; /**< @internal @brief image width */
int imHeight ; /**< @internal @brief image height */
int stepX ; /**< frame sampling step X */
int stepY ; /**< frame sampling step Y */
int boundMinX ; /**< frame bounding box min X */
int boundMinY ; /**< frame bounding box min Y */
int boundMaxX ; /**< frame bounding box max X */
int boundMaxY ; /**< frame bounding box max Y */
/** descriptor parameters */
VlDsiftDescriptorGeometry geom ;
int useFlatWindow ; /**< flag: whether to approximate the Gaussian window with a flat one */
double windowSize ; /**< size of the Gaussian window */
int numFrames ; /**< number of sampled frames */
int descrSize ; /**< size of a descriptor */
VlDsiftKeypoint *frames ; /**< frame buffer */
float *descrs ; /**< descriptor buffer */
int numBinAlloc ; /**< buffer allocated: descriptor size */
int numFrameAlloc ; /**< buffer allocated: number of frames */
int numGradAlloc ; /**< buffer allocated: number of orientations */
float **grads ; /**< gradient buffer */
float *convTmp1 ; /**< temporary buffer */
float *convTmp2 ; /**< temporary buffer */
} VlDsiftFilter ;主要逻辑代码注释如下：
vl_dsift.c中的mex函数：
void
mexFunction(int nout, mxArray *out[],
int nin, const mxArray *in[])
{
/*前面是输入参数的处理，此处从Do job部分开始*/
//参数初始化
int numFrames ;
int descrSize ;
VlDsiftKeypoint const *frames ;
float const *descrs ;
int k, i ;
VlDsiftFilter *dsift ;
//M,N是图像width,height
dsift = vl_dsift_new (M, N) ;
vl_dsift_set_geometry(dsift, &geom) ;
//step[0],step[1]分别代表x,y方向上的移动步长，verbose打开时值不相等
vl_dsift_set_steps(dsift, step[0], step[1]) ;
if (bounds) {
vl_dsift_set_bounds(dsift,
VL_MAX(bounds[1], 0),
VL_MAX(bounds[0], 0),
VL_MIN(bounds[3], M - 1),
VL_MIN(bounds[2], N - 1));
}
vl_dsift_set_flat_window(dsift, useFlatWindow) ;
if (windowSize >= 0) {
vl_dsift_set_window_size(dsift, windowSize) ;
}
numFrames = vl_dsift_get_keypoint_num (dsift) ;
descrSize = vl_dsift_get_descriptor_size (dsift) ;
geom = *vl_dsift_get_geometry (dsift) ;
/*此处省略一部分代码，是处理verbose打开的情况*/
//计算sift特征
vl_dsift_process (dsift, data) ;
//这里得到的frames中还包含norm
frames = vl_dsift_get_keypoints (dsift) ;
descrs = vl_dsift_get_descriptors (dsift) ;
/*后面将frames和descrs中的数据处理(归一化等)再输出*/
}vl_dsift_process函数在dsift.c中：
void vl_dsift_process (VlDsiftFilter* self, float const* im)
{
int t, x, y ;
_vl_dsift_alloc_buffers (self) ;
for (t = 0 ; t < self->geom.numBinT ; ++t)
memset (self->grads[t], 0,
sizeof(float) * self->imWidth * self->imHeight) ;
#undef at
#define at(x,y) (im[(y)*self->imWidth+(x)])
//对每一个像素点计算梯度(幅值和幅角)，norm，
for (y = 0 ; y < self->imHeight ; ++ y) {
for (x = 0 ; x < self->imWidth ; ++ x) {
float gx, gy ;
float angle, mod, nt, rbint ;
int bint ;
//y方向梯度
if (y == 0) {
gy = at(x,y+1) - at(x,y) ;
} else if (y == self->imHeight - 1) {
gy = at(x,y) - at(x,y-1) ;
} else {
gy = 0.5F * (at(x,y+1) - at(x,y-1)) ;
}
//x方向梯度
if (x == 0) {
gx = at(x+1,y) - at(x,y) ;
} else if (x == self->imWidth - 1) {
gx = at(x,y) - at(x-1,y) ;
} else {
gx = 0.5F * (at(x+1,y) - at(x-1,y)) ;
}
//计算幅角
angle = vl_fast_atan2_f (gy,gx) ;
//计算幅值
mod = vl_fast_sqrt_f (gx*gx + gy*gy) ;
//计算8个方向的值，把角度值转换成实数值
nt = vl_mod_2pi_f (angle) * (self->geom.numBinT / (2*VL_PI)) ;
bint = (int) vl_floor_f (nt) ;
rbint = nt - bint ;
//存梯度信息，统计直方图
self->grads [(bint ) % self->geom.numBinT][x + y * self->imWidth] = (1 - rbint) * mod ;
self->grads [(bint + 1) % self->geom.numBinT][x + y * self->imWidth] = ( rbint) * mod ;
}
}
//这里的flat_window是一种比高斯函数较快的平滑方法
if (self->useFlatWindow) {
_vl_dsift_with_flat_window(self) ;
} else {
_vl_dsift_with_gaussian_window(self) ;
}
{
VlDsiftKeypoint* frameIter = self->frames ;
float * descrIter = self->descrs ;
int framex, framey, bint ;
int frameSizeX = self->geom.binSizeX * (self->geom.numBinX - 1) + 1 ;
int frameSizeY = self->geom.binSizeY * (self->geom.numBinY - 1) + 1 ;
int descrSize = vl_dsift_get_descriptor_size (self) ;
float deltaCenterX = 0.5F * self->geom.binSizeX * (self->geom.numBinX - 1) ;
float deltaCenterY = 0.5F * self->geom.binSizeY * (self->geom.numBinY - 1) ;
float normConstant = frameSizeX * frameSizeY ;
for (framey = self->boundMinY ;
framey <= self->boundMaxY - frameSizeY + 1 ;
framey += self->stepY) {
for (framex = self->boundMinX ;
framex <= self->boundMaxX - frameSizeX + 1 ;
framex += self->stepX) {
frameIter->x = framex + deltaCenterX ;
frameIter->y = framey + deltaCenterY ;
//norm是以当前像素点为中心点的patch中所有像素的梯度幅值的平均值
{
float mass = 0 ;
for (bint = 0 ; bint < descrSize ; ++ bint)
mass += descrIter[bint] ;
mass /= normConstant ;
frameIter->norm = mass ;
}
/* L2 normalize */
_vl_dsift_normalize_histogram (descrIter, descrIter + descrSize) ;
/* clamp */
for(bint = 0 ; bint < descrSize ; ++ bint)
if (descrIter[bint] > 0.2F) descrIter[bint] = 0.2F ;
/* L2 normalize */
_vl_dsift_normalize_histogram (descrIter, descrIter + descrSize) ;
frameIter ++ ;
descrIter += descrSize ;
} /* for framex */
} /* for framey */
}
}时间： 2024-11-08 21:40:33