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 image
FastHessian::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 setsaveParameters(octaves, intervals, init_sample, thresh);
}//-------------------------------------------------------
//! 包含积分图的fasthessian结构构建
//! Constructor with image
FastHessian::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 setsaveParameters(octaves, intervals, init_sample, thresh);// Set the current imagesetIntImage(img);
}//-------------------------------------------------------FastHessian::~FastHessian()
{for (unsigned int i = 0; i < responseMap.size(); ++i){delete responseMap[i];}
}//-------------------------------------------------------
//！ 初始化参数
//! Save the parameters
void FastHessian::saveParameters(const int octaves, const int intervals, const int init_sample, const float thresh)
{// Initialise variables with bounds-checked valuesthis->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 imagethis->img = img;i_height = img->height;i_width = img->width;
}//-------------------------------------------------------//! 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 iptsipts.clear();// 建立尺度空间// Build the response mapbuildResponseMap();// 得到极值点// Get the response layersResponseLayer *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 spacefor (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 responses
void 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 layersfor(unsigned int i = 0; i < responseMap.size(); ++i)  delete responseMap[i];responseMap.clear();// 得到图像的参数// Get image attributesint w = (i_width / init_sample);//宽=原始图像宽/原始抽样倍数int h = (i_height / init_sample);//高=原始图像高/原始抽样倍数int s = (init_sample);//原始抽样倍数// 创建尺度空间所有层// Calculate approximated determinant of hessian valuesif (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 imagefor (unsigned int i = 0; i < responseMap.size(); ++i){buildResponseLayer(responseMap[i]);}
}//-------------------------------------------------------
//! 计算尺度空间每层的hessian值及laplacian值（及迹的值
//! Calculate DoH responses for supplied layer
void 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 coordinatesr = ar * step;c = ac * step; // 计算hessian成员值// Compute response componentsDxx = 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 sizeDxx *= inverse_area;Dyy *= inverse_area;Dxy *= inverse_area;// 保存// Get the determinant of hessian response & laplacian signresponses[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 responserl->coords.push_back(std::make_pair<int,int>(r,c));
#endif}}
}//-------------------------------------------------------
//! 极值点检测，r，c像素点坐标，t，m，b，尺度空间中连续的三层
//! Non Maximal Suppression function
int FastHessian::isExtremum(int r, int c, ResponseLayer *t, ResponseLayer *m, ResponseLayer *b)
{//计算边界，检测r/c是否越界// bounds checkint 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 maximumif (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 bottomint 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 extremumdouble xi = 0, xr = 0, xc = 0;interpolateStep(r, c, t, m, b, &xi, &xr, &xc );//插值法逼近真实极值点// 该点是否有效逼近真实极值点// If point is sufficiently close to the actual extremumif( 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;
}//-------------------------------------------------------