SURF算法:OpenSURF部分代码,分析与注释

OpenSURF里使用的默认参数(octaves=5,
intervals=4, init_sample=2, thresh=0.0004)表示生成的金字塔pyramid是5叠,每叠4层的,采样步长为2,在判断是否为极值(is
extrema)时阈值为0.0004.

OpenSURF中的主要文件有
fasthessian.h&fasthessian.cpp,integral.h&integral.cpp, ipoint.h&ipoint.cpp, Kmeans.h, main.cpp,responselayer.h,surf.cpp&surf.h, surflib.h和utils.cpp.下面主要是跟SURF论文有关的几个函数的代码和注释.

在surlib.h中可以看出算法的计算流程是:1.Integral求积分图;
2.FastHessian计算金字塔pyramid;
3.getIpoints检测关键点key
point或者兴趣点interestpoint;
4.getOrientation对每一个ipoint计算参考方向;
5.getDescriptors得到描述符.

ipoint.h:定义特征点

/***********************************************************

* --- OpenSURF --- *

* C. Evans, Research Into Robust Visual Features, *

* MSc University of Bristol,2008.
*************************************************************/

classIpoint {

//重载运算符计算两个描述符之间的欧几里得距离

floatoperator-(const Ipoint &rhs)

{

floatsum=0.f;

for(inti=0; i < 64; ++i)

sum +=(this->descriptor[i] - rhs.descriptor[i])*(this->descriptor[i] -rhs.descriptor[i]);

returnsqrt(sum);

};

//!Coordinates of the detected interest point兴趣点坐标

float x, y;

//!Detected scale兴趣点所在的尺度,在interpolateExtremum函数中通过插值求得.

floatscale;

//!Orientation measured anti-clockwise from +ve x-axis兴趣点的参考方向,从x轴正方向开始逆时针0到360度,算法里采用0到2pi.

floatorientation;

//! Sign oflaplacian for fast matching purposes计算fastheassian结果的拉普拉斯符号.

intlaplacian;

//! Vectorof descriptor components 64位描述符

float descriptor[64];

//!Placeholds for point motion (can be used for frame to frame motion analysis)两幅相片一对匹配点间坐标的差值,用于计算行对定位,摄站坐标系的移动.

float dx,dy;

//! Used tostore cluster index再度一张图像中特征点聚类分析后得到的聚类的编号.

intclusterIndex;

};

ipoint.cpp

//! Populate IpPairVec with matchedipts 实现寻找描述符距离最近的两个兴趣点

void getMatches(IpVec &ipts1, IpVec&ipts2, IpPairVec &matches)

{

float dist, d1, d2;

Ipoint *match;

matches.clear();

for(unsigned int i = 0; i < ipts1.size(); i++)

{

d1 = d2 = FLT_MAX;

for(unsigned int j = 0; j < ipts2.size(); j++)

{

dist = ipts1[i] - ipts2[j];

if(dist<d1) // if this feature matches better than current best

{//循环寻找更接近匹配点

d2 = d1;

d1 = dist;

match = &ipts2[j];

}

else if(dist<d2) // this feature matches better than second best

{

d2 = dist;

}

}

// If match has a d1:d2 ratio < 0.65 ipoints are a match跟SIFT一样,当最近的描述符远小于第二进的描述符时认为找到匹配点,个人感觉不如用KNN.

if(d1/d2 < 0.65)

{

// Store the change in position

ipts1[i].dx = match->x - ipts1[i].x;

ipts1[i].dy = match->y - ipts1[i].y;

matches.push_back(std::make_pair(ipts1[i], *match));

}

}

还有一部分是计算相对定位的代码,没有列出.

}

integral.h定义积分图

inline float BoxIntegral(IplImage *img,int row, int col, int rows, int cols)

{//求矩形ABCD中的积分=A+D-B-C;

float *data = (float *) img->imageData;

int step = img->widthStep/sizeof(float);

// The subtraction by one for row/col is because row/col is inclusive.

int r1 = std::min(row, img->height) - 1;

int c1 = std::min(col, img->width) - 1;

int r2 = std::min(row + rows, img->height) - 1;

int c2 = std::min(col + cols, img->width) - 1;

float A(0.0f), B(0.0f), C(0.0f), D(0.0f);

if (r1 >= 0 && c1 >= 0) A = data[r1 * step + c1];//坐标转换

if (r1 >= 0 && c2 >= 0) B = data[r1 * step + c2];

if (r2 >= 0 && c1 >= 0) C = data[r2 * step + c1];

if (r2 >= 0 && c2 >= 0) D = data[r2 * step + c2];

return std::max(0.f, A - B - C + D);

}

integral.cpp

//! Computes the integral image of imageimg. Assumes source image to be a对图形求积分图

IplImage *Integral(IplImage *source)

{//在这里作者求出的积分图,保存在了一个一维的数组里

// first row only//先计算第一行

float rs = 0.0f;

for(int j=0; j<width; j++)

{

rs += data[j];

i_data[j] = rs;

}

// remaining cells are sum above and to the left//积分图每一个点的值是左上部分的积分

for(int i=1; i<height; ++i)

{

rs = 0.0f;

for(int j=0; j<width; ++j)

{

rs += data[i*step+j];

i_data[i*step+j] = rs + i_data[(i-1)*step+j];

}

}

}

responselayer.h:对金字塔pyramid中的每一叠octave中每一层interval的定义

class ResponseLayer

{

public:

int width, height, step, filter;

float *responses;

unsigned char *laplacian;

ResponseLayer(int width, int height, int step, int filter)

{

assert(width > 0 && height > 0);//assert计算括号内条件,成立则继续执行,不成立则报告错误并终止运行

this->width = width;//interval的宽

this->height = height;//interval的高

this->step = step;//采样的步长

this->filter = filter;//threshold计算fasthessian时的box的尺寸

responses = new float[width*height];//以下申请空间

laplacian = new unsigned char[width*height];

memset(responses,0,sizeof(float)*width*height);

memset(laplacian,0,sizeof(unsigned char)*width*height);

}

inline unsigned char getLaplacian(unsignedint row, unsigned int column)

{//读取第
row行第column列的laplacian的值

return laplacian[row * width + column];

}

inline unsigned char getLaplacian(unsigned int row, unsigned int column,ResponseLayer *src)

{//src是参考的interval,读取调用这个函数的interval,对应于src坐标系的row,column处的
laplacian的值.

int scale = this->width / src->width;

#ifdef RL_DEBUG

assert(src->getCoords(row, column) == this->getCoords(scale * row,scale * column));

#endif

return laplacian[(scale * row) * width + (scale * column)];

}

另有 两个getResponse函数
与以上部分相同没有给出.

inline std::pair<int,int>getCoords(unsigned int row, unsigned int column)

{//返回interval中的第row行column列对应积分图的坐标,因为作者把积分图保存为了一个一维矩阵

return coords[row * width + column];

}

inline std::pair<int,int> getCoords(unsigned int row, unsigned intcolumn, ResponseLayer *src)

{//计算调用此函数的interval和被调用的interval之间的坐标转换.然后调用gerCoords

int scale = this->width / src->width;

return coords[(scale * row) * width + (scale * column)];

}

}

fasthessian.cpp

//! Find the image features and writeinto vector of feature

void FastHessian::getIpoints()

{//游历金字塔寻找兴趣点getIpoints和SIFT中的类似,不过
SIFT 中每一个octave中的层的大小是一样的.

//
filter index map//金字塔的组成,每叠包含的层,第一叠octave包含第(0,1,2,3)层...

staticconst int filter_map [OCTAVES][INTERVALS] = {{0,1,2,3}, {1,3,4,5}, {3,5,6,7},{5,7,8,9}, {7,9,10,11}};

// Clear the vector of exisiting ipts

ipts.clear();

// Build the response map

buildResponseMap();//建立金字塔,确定金字塔的结构

// Get the response layers

ResponseLayer *b, *m, *t;//实现金字塔,计算每一个具体值

for (int o = 0; o < octaves; ++o) for (int i = 0; i <= 1; ++i)

{//对每一个octave从小到上一次去三个层并对中间的每一个点层判定是否为极值点

b = responseMap.at(filter_map[o][i]);
//b: base底部

m = responseMap.at(filter_map[o][i+1]);//m: middle中间

t = responseMap.at(filter_map[o][i+2]);//t: top顶层

// loop over middle response layer at density of the most

// sparse layer (always top), to find maxima across scale and space

for (int r = 0; r < t->height; ++r)

{

for (int c = 0; c < t->width; ++c)

{

if (isExtremum(r, c, t, m, b))
//调用isExtremum函数与SIFT类似

{

interpolateExtremum(r, c, t, m, b);//调用interpolateExtremum函数与SIFT类似

}

}

}

}

}

//-------------------------------------------------------

//! Build map of DoH responses

void FastHessian::buildResponseMap()//定义金字塔的结构,对octave和interval申请

{

// Calculate responses for the first 4 octaves://对用不同的层,计算fasthessian时使用的不同boxsize

// Oct1: 9, 15, 21, 27

// Oct2: 15, 27, 39, 51

// Oct3: 27, 51, 75, 99

// Oct4: 51, 99, 147,195

// Oct5: 99, 195,291,387

// Get image attributes

int w = (i_width / init_sample);//宽

int h = (i_height / init_sample);//高

int s = (init_sample);//采样步长

// Calculate approximated determinant of hessian values

if (octaves >= 1)

{

responseMap.push_back(new ResponseLayer(w, h, s, 9));

responseMap.push_back(new ResponseLayer(w, h, s, 15));

responseMap.push_back(new ResponseLayer(w, h, s, 21));

responseMap.push_back(new ResponseLayer(w, h, s, 27));

}

if (octaves >= 2)

{

responseMap.push_back(new ResponseLayer(w/2, h/2, s*2, 39));

responseMap.push_back(new ResponseLayer(w/2, h/2, s*2, 51));

}

if (octaves >= 3)

{

responseMap.push_back(new ResponseLayer(w/4, h/4, s*4, 75));

responseMap.push_back(new ResponseLayer(w/4, h/4, s*4, 99));

}

if (octaves >= 4)

{

responseMap.push_back(new ResponseLayer(w/8, h/8, s*8, 147));

responseMap.push_back(new ResponseLayer(w/8, h/8, s*8, 195));

}

if (octaves >= 5)

{

responseMap.push_back(new ResponseLayer(w/16, h/16, s*16, 291));

responseMap.push_back(new ResponseLayer(w/16, h/16, s*16, 387));

}

}

//-------------------------------------------------------

//! Calculate DoH responses for suppliedlayer

voidFastHessian::buildResponseLayer(ResponseLayer *rl)//实现每一个层,对层中每一个位置计算fasthessian

{计算boxintegral部分需要的参数每个参数的意义可以看论文

float *responses = rl->responses; // response storage

unsigned char *laplacian = rl->laplacian; // laplacian sign storage

int step = rl->step; // step size for thisfilter

int b = (rl->filter - 1) / 2; // border for this filter

int l = rl->filter / 3; // lobe for this filter(filter size / 3)

int w = rl->filter; // filter size

float inverse_area = 1.f/(w*w); // normalisation factor

float Dxx, Dyy, Dxy;

for(int r, c, ar = 0, index = 0; ar < rl->height; ++ar)

{

for(int ac = 0; ac < rl->width; ++ac,index++)

{

// get the image coordinates//坐标转换

r = ar * step;

c = ac * step;

// Compute response components

Dxx = BoxIntegral(img, r - l + 1, c - b, 2*l - 1, w)

- BoxIntegral(img, r - l + 1, c - l / 2, 2*l - 1, l)*3;

Dyy = BoxIntegral(img, r - b, c - l + 1, w, 2*l - 1)

- BoxIntegral(img, r - l / 2, c - l + 1, l, 2*l - 1)*3;

Dxy = + BoxIntegral(img, r - l, c + 1, l, l)

+ BoxIntegral(img, r + 1, c - l, l,l)

- BoxIntegral(img, r - l, c - l, l,l)

- BoxIntegral(img, r + 1, c + 1, l,l);

// Normalise the filter responses with respect to their size

Dxx *= inverse_area;//归一化

Dyy *= inverse_area;

Dxy *= inverse_area;

// Get the determinant of hessian response & laplacian sign

responses[index] = (Dxx * Dyy - 0.81f * Dxy * Dxy);//论文里fasthessian的计算公式

laplacian[index] = (Dxx + Dyy >= 0 ? 1 : 0);

// create list of the image coords for each response

rl->coords.push_back(std::make_pair<int,int>(r,c));

}

}

}

//-------------------------------------------------------

int FastHessian::isExtremum(int r, intc, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b)

{

// bounds check

int layerBorder = (t->filter + 1) / (2 * t->step);

if (r <= layerBorder || r >= t->height - layerBorder || c <=layerBorder || c >= t->width - layerBorder)

return 0;//是否太靠近图像边界

// check the candidate point in the middle layer is above thresh

float candidate = m->getResponse(r, c, t);//取中间层m的对应于顶层t坐标(r,c)处的值

if (candidate < thresh)

return 0;

for (int rr = -1; rr <=1; ++rr)

{

for (int cc = -1; cc <=1; ++cc)

{

// if any response in 3x3x3 is greater candidate not maximum

if (//在3*3*3的邻域范围内判断它是否是极大值

t->getResponse(r+rr, c+cc) >= candidate ||

((rr != 0 || cc != 0) && m->getResponse(r+rr, c+cc, t) >=candidate) ||

b->getResponse(r+rr, c+cc, t) >= candidate

)

return 0;

}

}

return 1;

}

//-------------------------------------------------------

//! Interpolate scale-space extrema tosubpixel accuracy to form an image feature.

voidFastHessian::interpolateExtremum(int r, int c, ResponseLayer *t, ResponseLayer*m, ResponseLayer *b)

{//插值计算准确的极值,计算与SIFT相同,可以看Lowe的SIFT论文

// get the step distance between filters

// check the middle filter is mid way between top and bottom

int filterStep = (m->filter - b->filter);

assert(filterStep > 0 && t->filter - m->filter ==m->filter - b->filter);

// Get the offsets to the actual location of the extremum

double xi = 0, xr = 0, xc = 0;

interpolateStep(r, c, t, m, b, &xi,& xr, &xc );//调用差值函数,不像SIFT用了循环,只计算了一次

// If point is sufficiently close to the actual extremum

if( fabs( xi ) < 0.5f && fabs( xr ) < 0.5f && fabs( xc ) < 0.5f )

{//如果,便宜没有超出0.5,的范围

Ipoint ipt;

ipt.x = static_cast<float>((c + xc) * t->step);//兴趣点的x坐标

ipt.y = static_cast<float>((r + xr) * t->step);//兴趣点的y坐标

ipt.scale = static_cast<float>((0.1333f) * (m->filter + xi *filterStep));//兴趣点的尺度

ipt.laplacian = static_cast<int>(m->getLaplacian(r,c,t));

ipts.push_back(ipt);

}

}

//-------------------------------------------------------

//! Performs one step of extremuminterpolation.

void FastHessian::interpolateStep(int r,int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b,

double* xi,double* xr, double* xc )

{

CvMat* dD, * H, * H_inv, X;

double x[3] = { 0 };

dD = deriv3D( r, c, t, m, b );

H = hessian3D( r, c, t, m, b );

H_inv= cvCreateMat( 3, 3, CV_64FC1 );

cvInvert( H, H_inv, CV_SVD );//求H的逆矩阵,伪逆矩阵

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];

}

//-------------------------------------------------------

//! Computes the partial derivatives inx, y, and scale of a pixel.

CvMat* FastHessian::deriv3D(int r, intc, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b)

{//求三维偏导数

CvMat* dI;

double dx, dy, ds;

dx = (m->getResponse(r, c + 1, t) - m->getResponse(r, c - 1, t)) /2.0;

dy = (m->getResponse(r + 1, c, t) - m->getResponse(r - 1, c, t)) /2.0;

ds = (t->getResponse(r, c) - b->getResponse(r, c, t)) / 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;

}

//-------------------------------------------------------

//! Computes the 3D Hessian matrix for apixel.

CvMat* FastHessian::hessian3D(int r, intc, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b)

{//三维hessian矩阵

CvMat* H;

double v, dxx, dyy, dss, dxy, dxs, dys;

v = m->getResponse(r, c, t);

dxx = m->getResponse(r, c + 1, t) + m->getResponse(r, c - 1, t) -2 * v;

dyy = m->getResponse(r + 1, c, t) + m->getResponse(r - 1, c, t) -2 * v;

dss = t->getResponse(r, c) + b->getResponse(r, c, t) - 2 * v;

dxy = ( m->getResponse(r + 1, c + 1, t) - m->getResponse(r + 1, c- 1, t) -

m->getResponse(r - 1, c + 1, t) + m->getResponse(r - 1, c - 1, t)) / 4.0;

dxs = ( t->getResponse(r, c + 1) - t->getResponse(r, c - 1) -

b->getResponse(r, c + 1, t) + b->getResponse(r, c - 1, t) ) / 4.0;

dys = ( t->getResponse(r + 1, c) - t->getResponse(r - 1, c) -

b->getResponse(r + 1, c, t) + b->getResponse(r - 1, c, t) ) / 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;

}

//到这里就完成了对interestpoint的提取;以下进行计算descriptor

surf.cpp

const float pi = 3.14159f;

//! lookup table for 2d gaussian (sigma= 2.5) where (0,0) is top left and (6,6) is bottom right

const double gauss25 [7][7] = {//15*15的高斯矩阵的右下部分,其他的地方是对称的,但是好像和用matlab计算出来的有一点差别

0.02546481, 0.02350698, 0.01849125, 0.01239505, 0.00708017, 0.00344629,
0.00142946,

0.02350698, 0.02169968, 0.01706957, 0.01144208, 0.00653582, 0.00318132,
0.00131956,

0.01849125, 0.01706957, 0.01342740, 0.00900066, 0.00514126, 0.00250252,
0.00103800,

0.01239505, 0.01144208, 0.00900066, 0.00603332, 0.00344629, 0.00167749,
0.00069579,

0.00708017, 0.00653582, 0.00514126, 0.00344629, 0.00196855, 0.00095820,
0.00039744,

0.00344629, 0.00318132, 0.00250252, 0.00167749, 0.00095820, 0.00046640,
0.00019346,

0.00142946, 0.00131956, 0.00103800, 0.00069579, 0.00039744, 0.00019346,
0.00008024

};

//-------------------------------------------------------

//以下是几个计算重要调用的辅助型函数

/! Calculate the value of the 2dgaussian at x,y

inline float Surf::gaussian(int x, inty, float sig

{//计算二维的高斯权重

return (1.0f/(2.0f*pi*sig*sig)) * exp( -(x*x+y*y)/(2.0f*sig*sig));

}

//-------------------------------------------------------

//! Calculate the value of the 2dgaussian at x,y

inline float Surf::gaussian(float x,float y, float sig)

{

return 1.0f/(2.0f*pi*sig*sig) * exp( -(x*x+y*y)/(2.0f*sig*sig));

}

//------------------------------------------------------

//! Calculate Haar wavelet responses inx direction

inline float Surf::haarX(int row, intcolumn, int s)

{//计算row
column处小波基长度为s的x方向的小波变换,这里类似于求该点处s范围内的灰度变化

returnBoxIntegral(img, row-s/2, column, s, s/2)

-1 * BoxIntegral(img, row-s/2, column-s/2, s, s/2);

}

//-------------------------------------------------------

//! Calculate Haar wavelet responses iny direction

inline float Surf::haarY(int row, int column,int s)

{//计算row
column处小波基长度为s的y方向的小波变换

return BoxIntegral(img, row, column-s/2, s/2, s)

-1 * BoxIntegral(img, row-s/2, column-s/2, s/2, s);

}

//-------------------------------------------------------

//! Get the angle from the +ve x-axis ofthe vector given by (X Y)

float Surf::getAngle(float X, float Y)

{//求角度,单位为弧度,以x轴正方向为0,逆时针,0到2pi

if(X > 0 && Y >= 0)

return atan(Y/X);

if(X < 0 && Y >= 0)

return pi - atan(-Y/X);

if(X < 0 && Y < 0)

return pi + atan(Y/X);

if(X > 0 && Y < 0)

return 2*pi - atan(-Y/X);

return 0;

}

//---------------------------------------------------------------

//! Assign the supplied Ipoint anorientation

void Surf::getOrientation()

{

Ipoint *ipt = &ipts[index];

float gauss = 0.f, scale = ipt->scale;

const int s = fRound(scale), r = fRound(ipt->y), c =fRound(ipt->x);

std::vector<float> resX(109), resY(109), Ang(109);

const int id[] = {6,5,4,3,2,1,0,1,2,3,4,5,6};

int idx = 0;

// calculate haar responses for points within radius of 6*scale

for(int i = -6; i <= 6; ++i)

{//计算离中心点6S范围内的点的小波变换,每一个点的小波变换形成一个二维向量

for(int j = -6; j <= 6; ++j)

{

if(i*i + j*j < 36)

{

gauss = static_cast<float>(gauss25[id[i+6]][id[j+6]]);

// could use abs() rather than id lookup,but this way is faster

resX[idx] = gauss * haarX(r+j*s, c+i*s, 4*s);//加权

resY[idx] = gauss * haarY(r+j*s, c+i*s, 4*s);

Ang[idx] = getAngle(resX[idx], resY[idx]);//计算角度,调用getAngle函数

++idx;

}

}

}

// calculate the dominant direction

float sumX=0.f, sumY=0.f;

float max=0.f, orientation = 0.f;

float ang1=0.f, ang2=0.f;

// loop slides pi/3 window around feature point

for(ang1 = 0; ang1 < 2*pi; ang1+=0.15f)

{//检测主方向,即在一个方向附近的角度范围内,上一步求出的向量在这一个角度范围的和达到极值,者认为这个范围的中心为主方向

ang2 = ( ang1+pi/3.0f > 2*pi ? ang1-5.0f*pi/3.0f : ang1+pi/3.0f);

sumX = sumY = 0.f;

for(unsigned int k = 0; k < Ang.size(); ++k)

{//游历以上求出的向量

// get angle from the x-axis of the sample point

const float & ang = Ang[k];

// determine whether the point is within the window//每一个向量是否落在这个角度范围内

if (ang1 < ang2 && ang1 < ang && ang < ang2)//将这个范围内的向量求和

{

sumX+=resX[k];

sumY+=resY[k];

}

else if (ang2 < ang1 && ((ang > 0 && ang <ang2) || (ang > ang1 && ang < 2*pi) ))

{

sumX+=resX[k];

sumY+=resY[k];

}

}

// if the vector produced from this window is longer than all

// previous vectors then this forms the new dominant direction

if (sumX*sumX + sumY*sumY > max)
//是否达到极值

{

// store largest orientation

max = sumX*sumX + sumY*sumY;

orientation = getAngle(sumX, sumY);//计算主方向

}

}

//assign orientation of the dominant response vector

ipt->orientation = orientation;

}

//-------------------------------------------------------

//! Get the modified descriptor. SeeAgrawal ECCV 08

//! Modified descriptor contributed byPablo Fernandez

void Surf::getDescriptor(bool bUpright)

{

int y, x, sample_x, sample_y, count=0;

int i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;

float scale, *desc, dx, dy, mdx, mdy, co, si;

float gauss_s1 = 0.f, gauss_s2 = 0.f;

float rx = 0.f, ry = 0.f, rrx = 0.f, rry = 0.f, len = 0.f;

float cx = -0.5f, cy = 0.f; //Subregion centers for the 4x4 gaussianweighting

Ipoint *ipt = &ipts[index];

scale = ipt->scale;

x = fRound(ipt->x);

y = fRound(ipt->y);

desc = ipt->descriptor;

if (bUpright)// UPBRIGHT-SURF不考虑主方向,没有以主方向进行旋转

{

co = 1;

si = 0;

}

else

{//主方向在原坐标系的投影

co = cos(ipt->orientation);

si = sin(ipt->orientation);

}

i = -8;

//Calculate descriptor for this interest point

while(i < 12)

{

j = -8;

i = i-4;

cx += 1.f;

cy = -0.5f;

while(j < 12)

{

dx=dy=mdx=mdy=0.f;

cy += 1.f;

j = j - 4;

ix = i + 5;

jx = j + 5;

//旋转矩阵R=|cos(sita)
sin(sita)| R的逆矩阵inv(R)=RtR的转置.

// |-sin(sita) cos(sita)|

xs = fRound(x + ( -jx*scale*si + ix*scale*co));

ys = fRound(y + ( jx*scale*co + ix*scale*si));

for (int k = i; k < i + 9; ++k)

{

for (int l = j; l < j + 9; ++l)

{

//Get coords of sample point on the rotated axis

sample_x = fRound(x + (-l*scale*si + k*scale*co));

sample_y = fRound(y + ( l*scale*co + k*scale*si));

//Get the gaussian weighted x and y responses

gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5f*scale);

rx = haarX(sample_y, sample_x, 2*fRound(scale));//计算采样点的XY方向的

ry = haarY(sample_y, sample_x, 2*fRound(scale));//哈尔小波变换

//Get the gaussian weighted x and y responses on rotated axis

rrx = gauss_s1*(-rx*si + ry*co);//加权

rry = gauss_s1*(rx*co + ry*si);

dx += rrx;//四个统计特征:
dx

dy += rry;// dy

mdx += fabs(rrx);//abs(dx)

mdy += fabs(rry);//abs(dy)

}

}

//Add the values to the descriptor vector

gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f);

desc[count++] = dx*gauss_s2;//加权

desc[count++] = dy*gauss_s2;

desc[count++] = mdx*gauss_s2;

desc[count++] = mdy*gauss_s2;

len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy) * gauss_s2*gauss_s2;

j += 9;

}

i += 9;

}

//Convert to Unit Vector

len = sqrt(len);

for(int i = 0; i < 64; ++i)

desc[i] /= len; //归一化

}

//----------------------------------------------------

到这里特征描述符也计算完了.在网上没有找到,SURF算法代码实现的说明.就去看了一下OpenSURF的代码,里面有很多与SIFT相似之处.

能力有限,对代码和算法的理解还有很多不足和不正确的地方,希望大家帮忙指正指教,也希望我的工作对您理解SURF过程有帮助.