surf源码解析，附加源码下载链接

说是surf，其实原算法是基于Notes on the OpenSURF Library这篇文档。

 

下载地址：http://opensurf1.googlecode.com/files/OpenSURFcpp.zip

 

首先定义 #define PROCEDURE 1

紧接着进入main()主函数：

?

int main(void)

{

  if(PROCEDURE == 1) returnmainImage();//静态图像的特征点监测

  if(PROCEDURE == 2) returnmainVideo();//通过摄像头实时监测，提取特征点

  if(PROCEDURE == 3) returnmainMatch();//图像匹配算法，在视频序列里寻找一个已知物体图像

  if(PROCEDURE == 4) returnmainMotionPoints();//视频中行为监测

  if(PROCEDURE == 5) returnmainStaticMatch();//静态图像间的特征点匹配

  if(PROCEDURE == 6) returnmainKmeans();//静态图像的k-means聚类

}

 我们以监测静态图像特征点（即1）为例了解surf算法如何提取特征点。

?

int mainImage(void)

{

  // Declare Ipoints and other stuff

  IpVec ipts;

 

  IplImage *img=cvLoadImage("imgs/sf.jpg");

 

  // Detect and describe interest points in the image

  clock_t start = clock();

  surfDetDes(img, ipts,false, 5, 4, 2, 0.0004f);//URF描述子特征提取实现函数

  clock_t end = clock();

 

  std::cout<<"OpenSURF found: " << ipts.size() << " interest points"<< std::endl;

  std::cout<<"OpenSURF took: " << float(end - start) / CLOCKS_PER_SEC  << " seconds"<< std::endl;

 

  // Draw the detected points

  drawIpoints(img, ipts);

   

  // Display the result

  showImage(img);

 

  return0;

}

EXP:

IpVec的定义在ipoint.h里，typedef std::vector<Ipoint> IpVec; 也就是说，IpVec是一个vector向量，每个成员是一个Ipoint。而Ipoint是一个类，也在ipoint.h里，作者将它按照结构体使用，都是public。

clock()的使用是为了测试程序运行的时间，一个记录起始时间，一个记录结束时间，最终将两者只差输出，即surfDetDes()特征提取所需要的时间。

drawIpoints()函数是在img基础上增加特征点的描述，说白了，就是添加特征点的函数。

 

那么，现在进入到surfDetDes()函数内部来探个究竟吧。

surfDetDes定义在surflib.h里面：

?

//! Library function builds vector of described interest points

inline void surfDetDes(IplImage *img, /* image to find Ipoints in */

                       std::vector<Ipoint> &ipts,/* reference to vector of Ipoints */

                       bool upright =false,/* run in rotation invariant mode? */

                       int octaves = OCTAVES,/* number of octaves to calculate */

                       int intervals = INTERVALS,/* number of intervals per octave */

                       int init_sample = INIT_SAMPLE,/* initial sampling step */

                       float thres = THRES/* blob response threshold */)

{

  // Create integral-image representation of the image

  IplImage *int_img = Integral(img);

   

  // Create Fast Hessian Object

  FastHessian fh(int_img, ipts, octaves, intervals, init_sample, thres);

  

  // Extract interest points and store in vector ipts

  fh.getIpoints();

   

  // Create Surf Descriptor Object

  Surf des(int_img, ipts);

 

  // Extract the descriptors for the ipts

  des.getDescriptors(upright);

 

  // Deallocate the integral image

  cvReleaseImage(&int_img);

}

　　总的来说，surfDetDes()里面规划了具体的步骤：

1.获取积分图像init_img;

2.计算Hessian矩阵特征fh;

3.利用fh提取特征点;

4.创建surf特征描述符，并和特征点贯联;

5.释放积分图像。

 

 

---

①积分图像的实现

 首先在Integral.cpp里面找到Integral()，如下：

View Code?

IplImage *Integral(IplImage *source)

{

  // convert the image to single channel 32f

  IplImage *img = getGray(source);

  IplImage *int_img = cvCreateImage(cvGetSize(img), IPL_DEPTH_32F, 1);

 

  // set up variables for data access

  int height = img->height;

  int width = img->width;

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

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

  float *i_data = (float *) int_img->imageData; 

 

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

    }

  }

 

  // release the gray image

  cvReleaseImage(&img);

 

  // return the integral image

  returnint_img;

}

1.　　首先将原输入转化为灰度图像，并创建一个大小等于灰度图像gray-image的图像数组--积分图像int_img。

   2.　　获取图像的信息，比如大小（高height和宽width）以及gray-image和积分图像int_img的数据首地址data && i_data。(注意此时数据类型为float)

      3.　　首先计算第一行像素的积分值，相当于一维数据的累加。其他数据通过迭代计算获取，i_data[i*step+j] = rs + i_data[(i-1)*step+j]，若当前点为（i₀,j₀）,其中rs就为第（i+1）行（即i=i₀）行j<=j₀之前所有像素值和。 如下所示：

                                         [其中黑色为当前点i_data[i*step+j]，绿色为i_data[(i-1)*step+j]，而rs=横放着的和黑色点同行的那块矩形框对应的区域像素值之和]

             4.　　释放灰度图像，并返回积分图像。

 

相关函数integral.h中的BoxIntegral().

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

其中，几个参数意思分别为源图像，row,col为A点的坐标值，rows和cols分别为高和宽。

利用上面的积分图像计算   A      B   这样的box区域里面所有像素点的灰度值之和。S=int_img(D)+int_img(A)-int_img(B)-int_img(C).

                                  C      D 

---

②Hessian矩阵特征的计算

FastHessian，计算hessian矩阵的类，它的定义在fasthessian.h里，实现在fasthessian.cpp里。

View Code?

class FastHessian {

   

  public:

    

    //! Constructor without image

    FastHessian(std::vector<Ipoint> &ipts,

                const int octaves = OCTAVES,

                const int intervals = INTERVALS,

                const int init_sample = INIT_SAMPLE,

                const float thres = THRES);

 

    //! Constructor with image

    FastHessian(IplImage *img,

                std::vector<Ipoint> &ipts,

                const int octaves = OCTAVES,

                const int intervals = INTERVALS,

                const int init_sample = INIT_SAMPLE,

                const float thres = THRES);

 

    //! Destructor

    ~FastHessian();

 

    //! Save the parameters

    void saveParameters(const int octaves,

                        const int intervals,

                        const int init_sample,

                        const float thres);

 

    //! Set or re-set the integral image source

    void setIntImage(IplImage *img);

 

    //! Find the image features and write into vector of features

    void getIpoints();

     

  private:

 

    //---------------- Private Functions -----------------//

 

    //! Build map of DoH responses

    void buildResponseMap();

 

    //! Calculate DoH responses for supplied layer

    void buildResponseLayer(ResponseLayer *r);

 

    //! 3x3x3 Extrema test

    int isExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b);   

     

    //! Interpolation functions - adapted from Lowe's SIFT implementation

    void interpolateExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b);

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

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

    CvMat* deriv3D(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b);

    CvMat* hessian3D(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b);

 

    //---------------- Private Variables -----------------//

 

    //! Pointer to the integral Image, and its attributes

    IplImage *img;

    int i_width, i_height;

 

    //! Reference to vector of features passed from outside

    std::vector<Ipoint> &ipts;

 

    //! Response stack of determinant of hessian values

    std::vector<ResponseLayer *> responseMap;

 

    //! Number of Octaves

    int octaves;

 

    //! Number of Intervals per octave

    int intervals;

 

    //! Initial sampling step for Ipoint detection

    int init_sample;

 

    //! Threshold value for blob resonses

    float thresh;

};

　　在public里面定义了两种构造函数分别对应有无源图像这一项参数，紧接着还定义了析构函数~FastHessian等函数。下面在fasthessian.cpp对这些函数的实现一一解释：

两个构造函数都是调用了saveParameters(octaves, intervals, init_sample, thresh)设置构造金字塔的参数，而带图像的构造函数另外多加了一句setIntImage(img)用来设置当前图像。

 

?

//! Save the parameters

void FastHessian::saveParameters(const int octaves, const int intervals,

                                 const int init_sample, const float thresh)

{

  // Initialise variables with bounds-checked values

  this->octaves =

    (octaves > 0 && octaves <= 4 ? octaves : OCTAVES);

  this->intervals =

    (intervals > 0 && intervals <= 4 ? intervals : INTERVALS);

  this->init_sample =

    (init_sample > 0 && init_sample <= 6 ? init_sample : INIT_SAMPLE);

  this->thresh = (thresh >= 0 ? thresh : THRES);

}

 

//! Set or re-set the integral image source

void FastHessian::setIntImage(IplImage *img)

{

  // Change the source image

  this->img = img;

 

  i_height = img->height;

  i_width = img->width;

}

　　由于在h头文件中已设置

?

static const int OCTAVES = 5;//组数

static const int INTERVALS = 4;//每组层数

static const float THRES = 0.0004f;//阈值

static const int INIT_SAMPLE = 2;//初始采样因子

　 所以　saveParameters的作用就是调整参数，以防过大或过小。

FastHessian::getIpoints()提取兴趣点：

View Code?

//! Find the image features and write into vector of features

void FastHessian::getIpoints()

{

  // filter index map

  static const 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

  ...<br>  ...

}

　　首先初始化filter_map，清空标记特征点的ipts结构体。

     创建高斯平滑层函数参数ResponseMap()，大小与论文所给完全一致，

 　　　　　　　　// 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

这些都是每组模板的大小，每组间隔递增，6,12,24,48,96 。ResponseMap这个结构体向量包含4个参数ResponseLayer(int width, int height, int step, int filter)定义在responsemap.h里面，其中width和height等于实际图像大小除以step(step初始值为2)，而filter则是滤波器半径。

 

然后使用buildResponseLayer(responseMap[i])对图像处理后将数据存放在responses和laplacian两个数组里面。

View Code?

void FastHessian::buildResponseLayer(ResponseLayer *rl)

{

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

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

  int step = rl->step;                     // step size for this filter  滤波器尺度因子

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

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

 

#ifdef RL_DEBUG

      // create list of the image coords for each response

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

#endif

    }

  }

}

　　其中计算Dxy和Dyy的示意图如下，其他的Dxx(Dyy的转置)读者自己参考。【此时w=9,l=9/3=3,对应于右图的总计算区域高度和宽度2*l-1】

 

                                                                                                        

  圆点为当前点，将红色方形区域1内的积分值减去绿色方形2区域内的积分值=Dxy*w²                  绿色方形区域2内的积分值减去2*红色色方形区域1内的积分值=Dyy*w² (==用一整块区域-3*红色区域）

最后将计算后的结果存放到ResponseLayer里面的response和laplacian一维数组里，数组的大小即为ResponseLayer->width * ResponseLayer->width。

 

这样就计算出了每一层的所有像素点的det(Happrox)=Dxx*Dyy-(0.9*Dxy)²,下面开始判断当前点是否是极值点。

 

---

③根据上一步计算每层的每个点的det(Happrox)，判断极值点。

在fasthessian.cpp里面找到getIpoints()，下面开始抽取每组（共octaves组）相邻的3层（共有4层，每组有2个这样层的满足）：                            

View Code?

ResponseLayer *b, *m, *t;

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

  {

    b = responseMap.at(filter_map[o][i]);

    m = responseMap.at(filter_map[o][i+1]);

    t = responseMap.at(filter_map[o][i+2]);

 

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

    // sparse layer (always top), to find maxima across scale and space<br>    //总是在最上层去寻找极大值点

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

    {

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

      {

        if(isExtremum(r, c, t, m, b))

        {

          interpolateExtremum(r, c, t, m, b);

        }

      }

    }

  }

　　在判断的时候用到了isExtremum()和interpolateExtremum()子函数，在当前cpp内找到它们，并分析。

View Code?

//! Non Maximal Suppression function

int FastHessian::isExtremum(int r, int c, 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)

    return0;

 

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

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

  if(candidate < thresh) //小于某一阈值也认定为不是极值点

    return0;

 

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

  {

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

    {

      // if any response in 3x3x3 is greater candidate not maximum<br>      //只要3x3x3邻域存在一点大于当前点的response值，只可认定该点非极值点

      if(

        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

        )//优先级&&>||

        return0;

    }

  }

 

  return1;

}

　　大家务必注意，由于各层大小不一致，顾在比较大小的时候，要统一到统一尺度下，在getResponse()里面有所体现，即先计算两者尺度之比，比如尺度之比为2，最上层当前点为（20，25），那么需要比较的紧邻下层点为（40,50）而不是下层的（20,25）。再看若当前点为极值点，那么调用interpolateExtremum():

View Code?

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

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

{

  // 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 );

 

  // If point is sufficiently close to the actual extremum

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

  {

    Ipoint ipt;

    ipt.x = static_cast<float>((c + xc) * t->step);

    ipt.y = static_cast<float>((r + xr) * t->step);

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

  }

}

　　本函数实现功能，利用插值精确计算极值点在原图像中的位置，并保存在ipt中。

      打开interpolateStep，其中deriv3D是求当前点的3的方向上的一阶导，hessian3D是求当前点的3维hessian二阶导，最后计算出3个方向的偏差值xi,xr,xc .

View Code?

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 );//用svd法求逆矩阵H_inv

  cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );//新建X

  cvGEMM( H_inv, dD, -1, NULL, 0, &X, 0 );//X=-dD*H_inv

 

  cvReleaseMat( &dD );

  cvReleaseMat( &H );

  cvReleaseMat( &H_inv );

 

  *xi = x[2];

  *xr = x[1];

  *xc = x[0];

}

　　下面开始在特征点周围提取特征描述符

---

④ 提取特征描述符

转到surf.cpp里寻找getDescriptors()，由于upright初始化设置为false（为特征点分配主方向，并旋转找到邻域提取特征描述符），若为true,则直接提取邻域特征描述符。

View Code?

//! Describe all features in the supplied vector

void Surf::getDescriptors(bool upright)

{

  // Check there are Ipoints to be described

  if(!ipts.size()) return;

 

  // Get the size of the vector for fixed loop bounds

  int ipts_size = (int)ipts.size();

 

  if(upright)

  {

    // U-SURF loop just gets descriptors

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

    {

      // Set the Ipoint to be described

      index = i;

 

      // Extract upright (i.e. not rotation invariant) descriptors

      getDescriptor(true);

    }

  }

  else

  {

    // Main SURF-64 loop assigns orientations and gets descriptors

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

    {

      // Set the Ipoint to be described

      index = i;

 

      // Assign Orientations and extract rotation invariant descriptors

      getOrientation();

      getDescriptor(false);

    }

  }

}

　　其实两者区别在于提取特征描述符之前判断upright的时候是否需要多调用一句getOrientation()就是了。前者不考虑旋转，两幅图只是尺度有异，而后者还考虑了旋转不变性，更加全面。

 几个地方：   

1.如论文提出的采取r=6的圆域，共计包含109个像素点，大家可以算算，在此不多解释了。

2.像素点的方向由harrx和harry计算得到,相当于梯度公式里面的dx和dy，然后利用getAngle得到θ（分布在0~2∏，而不是简单的tan^-1(harrx/harry)取值在0~∏）。

?

//! Calculate Haar wavelet responses in x direction

inline float Surf::haarX(int row, int column, int 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 in y direction

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

{

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

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

}

 

float Surf::getAngle(float X, float Y)//计算每个点的方向角

{

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

    returnatan(Y/X);

 

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

    returnpi - atan(-Y/X);

 

  if(X < 0 && Y < 0)

    returnpi + atan(Y/X);

 

  if(X > 0 && Y < 0)

    return2*pi - atan(-Y/X);

 

  return0;

}

　　图示如右：                            （黑色代表-1，白色为1，即白色区域积分值减去黑色区域积分值）

 

在getOrienation()里面resx和resy分别是上面计算出来的harrx和harry分别乘上高斯模板系数gauss25。

?

void Surf::getOrientation()

{

  ......

   //将像素点按方向角分配到6个宽为60度的区间去，比如说可以将80度分配到ang1=90度，因为90-30<80<90+30

   // 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);//保证ang1+60 不会大于360度

    sumX = sumY = 0.f;

    for(unsigned int k = 0; k < Ang.size(); ++k)//对于半径为6的圆邻域，这里面的Ang.size()=109

    {

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

      }

      elseif (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)//寻找一个ang1使得角度在此区间的长度最大 ，然后计算出累加的resx和resy，得出的方向角

    {

      // store largest orientation

      max = sumX*sumX + sumY*sumY;

      orientation = getAngle(sumX, sumY);

    }

  }

  // assign orientation of the dominant response vector

  ipt->orientation = orientation;

}

　　好了，终于要开始提取特征描述符了哈~.~

View Code?

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 gaussian weighting

 

  Ipoint *ipt = &ipts[index];

  scale = ipt->scale;//尺度

  x = fRound(ipt->x);

  y = fRound(ipt->y);//空间位置 

  desc = ipt->descriptor;

 

  if(bUpright)

  {//不需旋转的情况

    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;//ix=i0+1,这里面i0和j0分别取得的值为-8，-3，2，7

      jx = j + 5;//iy=j0+1

 

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

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

      //fRound为4舍五入算法，最近邻插值寻找旋转对应点

 

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

          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;

          dy += rry;

          mdx += fabs(rrx);

          mdy += fabs(rry);

 

        }

      }

 

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

 

}

　　其中i和j分别取的值为-8，-3，2，7，很明显i,j确定的邻域为7-（-8）+1=16，16x16的邻域，旋转对应的在原图像的点位置为   

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

                               co = cos(ipt->orientation);                                                  si = sin(ipt->orientation);                        【ipt->orientation为特征点的方向角】

共有 4x4个子块（9x9=81个像素点），每个子块分别计算了其中16个dx,dy,|dx|,|dy|之和（当然还要考虑高斯滤波权重系数），则最终的特征描述符为4x4x4=64维向量。

 

main.cpp内mainImage函数内部drawIpoints(img, ipts)就不用再做解释了吧。

 

搞了一天，终于搞完了~~

 来个截图~~

 

转载自：blue_lghttp://www.cnblogs.com/blue-lg/archive/2012/07/20/2600436.html