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