opensurf具体过程分析

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

具体函数调用即surflib里面的surfDetDes函数。本文就如下代码段作个分析

 void surfDetDes(IplImage *img,  /* 检测兴趣点的图像 image to find Ipoints in */
					   std::vector<Ipoint> &ipts, /* 兴趣点向量的引用reference to vector of Ipoints */
					   bool upright, /*是否是旋转不变模式 run in rotation invariant mode? */
					   int octaves, /* 高斯金字塔的组number of octaves to calculate */
					   int intervals, /* 高斯金字塔每组的层number of intervals per octave层 */
					   int init_sample,  /* 初始抽样步骤 initial sampling step */
					   float 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();

	//创建Surf描述子对象  Create Surf Descriptor Object
	Surf des(int_img, ipts);

	// 提取ipts的描述子 Extract the descriptors for the ipts
	des.getDescriptors(upright);

	// Deallocate the integral image
	cvReleaseImage(&int_img);
}

首先是创建积分图像（IplImage *Integral(IplImage *source)），也就是说在积分图像中的每个像素点的值即该点的积分值，也就是以该点为右下角，图像原点为左上角的矩形区域内所有像素的总和。

积分过程具体如下，表1假设为 data数据，表2为需要生成的积分图i_data，首先对 i_data 第一行赋值，即相同位置 data同行 前面像素累加。对于第二行，则为相同位置data同行 前面像素累加再 加上 i_data中上一行同列的像素值。

表1

0	1	2	3
4	5	6	7

表2

0	0+1=1	0+1+2=3	0+1+2+3=6
4	5+4+1=10	6+5+4+3=18	7+6+5+4+6=28

积分图构造完毕后，接着是创建快速海森矩阵FastHessian
首先调用构造函数创建对象fh，主要是保证参数有效性以及初始化图像参数。

主要提取特征点函数：fh.getIpoints()，主要过程如下：

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下面判断是否是极值点
  ResponseLayer *b, *m, *t;
  for (int i = 0; i < octaves; ++i) //第几组
	  for (int j = 0; j <= 1; ++j) //第几列
	  {
		b = responseMap.at(filter_map[i][j]);//0指向filter9
		m = responseMap.at(filter_map[i][j+1]);//15
		t = responseMap.at(filter_map[i][j+2]);//21

		// loop over middle response layer at density of the most 
		// sparse layer (always top), to find maxima across scale and space
		for (int h = 0; h < t->height; ++h)
		{
		  for (int w = 0; w < t->width; ++w)
		  {
			if (isExtremum(h, w, t, m, b))
			{
			  interpolateExtremum(h, w, t, m, b);
			}
		  }
		}
	  }
}

首先初始化filter_map，如下所示共12个filtersize，与 filter_map[][]正好对应
9,15,21,27,39,51,75,99,147,195,291,387

然后调用buildResponseMap()：

主要执行操作responseMap.push_back(new ResponseLayer(w,   h,   s,   9));  共执行12次，即每种size的fiter。

ResponseLayer类的构造函数如下：

 ResponseLayer(int width, int height, int step, int filter)
  {
    assert(width > 0 && height > 0);
    
    this->width = width;
    this->height = height;
    this->step = step;//采样因子
    this->filter = filter;//盒子滤波的大小

    responses = new float[width*height];//存放每个大小的filer下的图像像素点的海森矩阵的行列式值
    laplacian = new unsigned char[width*height];//存放laplacian标识符，判断是暗背景上的亮斑还是亮背景的暗斑

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

通过buildResponseLayer(ResponseLayer * rl)函数计算海森矩阵的值

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
	  //整个区域-3倍中间区域，即权值为2的区域
      Dxx = BoxIntegral(img, r - l + 1, c - b, 2*l - 1, w)//BoxIntegral(积分图，起点行，起点列，高度，宽度)：求矩形区域内的积分，即像素总和A+D-B-C
          - 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);//Laplacian标识符
}

其中Dxx、Dyy、Dxy的计算说明如下，以filtersize=9为例，b=4，l=3，w=9，假设中心红点处的坐标值为为(r,c)=(4,4)，则整体区域的起点为(0,2)

Dyy=整个区域（2个黄色+绿色）-3*绿色

这样就计算出了每一层的所有像素点的det(Happrox)=dxx*dyy-(0.9*dxy)2，下面开始判断当前点是否是极值点。

代码部分回到getIpoints()4个for循环嵌套，抽取每组（共octaves组）相邻的3层(每组分别为4层)，所以总共能抽取2次，即i=0;i<=1的情况。

在判断极值点的时候分别用到了isExtremum()：非最大值抑制和interpolateExtremum()两个函数，