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 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 set saveParameters(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 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 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; } //------------------------------------------------------- //! 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 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 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 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 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 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 function int 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; } //-------------------------------------------------------