SURF源码分析之fasthessian.h和fasthessian.cpp

fasthessian.h和fasthessian.h主要完成关键点的检测

1、建立fasthesian结构

2、初始化fasthessian结构

3、建立DOH

4、找出极值点

5、确定关键点准确位置

源码分析：

/*********************************************************** *  --- OpenSURF ---                                       **  This library is distributed under the GNU GPL. Please   **  use the contact form at http://www.chrisevansdev.com    **  for more information.                                   **                                                          **  C. Evans, Research Into Robust Visual Features,         **  MSc University of Bristol, 2008.                        **                                                          *************************************************************/#include "integral.h"#include "ipoint.h"#include "utils.h"#include <vector>#include "responselayer.h"#include "fasthessian.h"using namespace std;//-------------------------------------------------------//! 不包含积分图的fasthessian结构构建//! Constructor without imageFastHessian::FastHessian(std::vector<Ipoint> &ipts,                          const int octaves, const int intervals, const int init_sample,                          const float thresh)                          : ipts(ipts), i_width(0), i_height(0){  // Save parameter set  saveParameters(octaves, intervals, init_sample, thresh);}//-------------------------------------------------------//! 包含积分图的fasthessian结构构建//! Constructor with imageFastHessian::FastHessian(IplImage *img, std::vector<Ipoint> &ipts,                          const int octaves, const int intervals, const int init_sample,                          const float thresh)                          : ipts(ipts), i_width(0), i_height(0){  // Save parameter set  saveParameters(octaves, intervals, init_sample, thresh);  // Set the current image  setIntImage(img);}//-------------------------------------------------------FastHessian::~FastHessian(){  for (unsigned int i = 0; i < responseMap.size(); ++i)  {    delete responseMap[i];  }}//-------------------------------------------------------//！ 初始化参数//! Save the parametersvoid 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 sourcevoid FastHessian::setIntImage(IplImage *img){  // Change the source image  this->img = img;  i_height = img->height;  i_width = img->width;}//-------------------------------------------------------//! Find the image features and write into vector of featuresvoid 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 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]);// 取连续的三层，计算最上层的极值点，将最上层的每个像素点与邻近的26个像素点比较    // 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))// 是否为极值点        {          interpolateExtremum(r, c, t, m, b);//用插值法计算精确的极值点位置        }      }    }  }}//-------------------------------------------------------//! 构建尺度空间//! Build map of DoH responsesvoid FastHessian::buildResponseMap(){  // 高斯滤波核前四组尺寸大小  // Calculate responses for the first 4 octaves:  // 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  // 清除已经存在的层  // Deallocate memory and clear any existing response layers  for(unsigned int i = 0; i < responseMap.size(); ++i)      delete responseMap[i];  responseMap.clear();  // 得到图像的参数  // 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));  }  // 提取每一层的hessian值及laplacian值  // Extract responses from the image  for (unsigned int i = 0; i < responseMap.size(); ++i)  {    buildResponseLayer(responseMap[i]);  }}//-------------------------------------------------------//! 计算尺度空间每层的hessian值及laplacian值（及迹的值//! Calculate DoH responses for supplied layervoid FastHessian::buildResponseLayer(ResponseLayer *rl){  float *responses = rl->responses;         // response storage hessian值存储数组  unsigned char *laplacian = rl->laplacian; // laplacian sign storage laplacian值存储数组  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;                      // hessian矩阵中四个元素，Dxy和Dyx值一样  //计算每个像素点的hessian值及laplacian值  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;   // 计算hessian成员值      // 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    }  }}  //-------------------------------------------------------//! 极值点检测，r，c像素点坐标，t，m，b，尺度空间中连续的三层//! Non Maximal Suppression functionint FastHessian::isExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b){  //计算边界，检测r/c是否越界  // 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;  // 检测hessian值是否大于阀值，如果小于，返回0  // check the candidate point in the middle layer is above thresh   float candidate = m->getResponse(r, c, t);  if (candidate < thresh)     return 0;   // 与附近26像素比较  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 (        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;}//-------------------------------------------------------//插值法逼近极值点准确位置，r、c像素点坐标，t，m，b尺度空间连续三层//! 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);  }}//-------------------------------------------------------//! 插值法算出极值点偏差//! Performs one step of extremum interpolation. 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 );//三维hessian矩阵计算  H_inv = cvCreateMat( 3, 3, CV_64FC1 );//创建64位单精度3*3矩阵  cvInvert( H, H_inv, CV_SVD );//转换矩阵格式  cvInitMatHeader( &X, 3, 1, CV_64FC1, x, CV_AUTOSTEP );//创建3*1矩阵头  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 in x, y, and scale of a pixel.CvMat* FastHessian::deriv3D(int r, int c, 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;}//-------------------------------------------------------//! 三维hessian矩阵计算//! Computes the 3D Hessian matrix for a pixel.CvMat* FastHessian::hessian3D(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b){  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;}//-------------------------------------------------------