[VLFeat]Dense Sift的C源码学习

VLFeat是一个很好用的开源库，其中实现了计算机视觉常用的算法，使用的语言是C和matlab。 官网： 在官网下载最新版本后，在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_{ s ; /**< scale */ double norm ; /**< SIFT descriptor norm */} VlDsiftKeypoint ;/** @brief Dense SIFT descriptor geometry */typedef struct VlDsiftDescriptorGeometry_{ numBinY ; binSizeY ; /**< size of bins along Y */} VlDsiftDescriptorGeometry ;typedef struct VlDsiftFilter_{ imHeight ;stepX ;boundMinX ;boundMaxX ;VlDsiftDescriptorGeometry geom ; numFrames ;/**< number of sampled frames */ int descrSize ;/**< size of a descriptor */ VlDsiftKeypoint *frames ; numBinAlloc ;numGradAlloc ;*convTmp1 ;/**< temporary buffer */ float *convTmp2 ;/**< temporary buffer */} VlDsiftFilter ;

主要逻辑代码注释如下： vl_dsift.c中的mex函数：

voidmexFunction(int nout, mxArray *out[],int nin, const mxArray *in[]){numFrames ;int descrSize ;VlDsiftKeypoint const *frames ;float const *descrs ;int k, i ; VlDsiftFilter *dsift ;//M,N是图像width,heightdsift = 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中还包含normframes = 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 */ }}

每个人在他的人生发轫之初，总有一段时光，