reference
一点补充
关于插值:利用了Hessian矩阵中不定点的特性,如果DOG空间某点为极值点,那么该点的Hessian存在特征点(不动点)。可以通过迭代进行插值计算。
static bool adjustLocalExtrema( const vector <Mat> & dog_pyr, KeyPoint& kpt, int octv,
int & layer, int & r, int & c, int nOctaveLayers,
float contrastThreshold, float edgeThreshold, float sigma )
{
const float img_scale = 1.f /(255 *SIFT_FIXPT_SCALE);
const float deriv_scale = img_scale*0.5f ;
const float second_deriv_scale = img_scale;
const float cross_deriv_scale = img_scale*0.25f ;
float xi=0 , xr=0 , xc=0 , contr=0 ;
int i = 0 ;
for ( ; i < SIFT_MAX_INTERP_STEPS; i++ )
{
int idx = octv*(nOctaveLayers+2 ) + layer;
const Mat& img = dog_pyr[idx];
const Mat& prev = dog_pyr[idx-1 ];
const Mat& next = dog_pyr[idx+1 ];
Vec3f dD((img.at<sift_wt>(r, c+1 ) - img.at<sift_wt>(r, c-1 ))*deriv_scale,
(img.at<sift_wt>(r+1 , c) - img.at<sift_wt>(r-1 , c))*deriv_scale,
(next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale);
float v2 = (float )img.at<sift_wt>(r, c)*2 ;
float dxx = (img.at<sift_wt>(r, c+1 ) + img.at<sift_wt>(r, c-1 ) - v2)*second_deriv_scale;
float dyy = (img.at<sift_wt>(r+1 , c) + img.at<sift_wt>(r-1 , c) - v2)*second_deriv_scale;
float dss = (next.at<sift_wt>(r, c) + prev.at<sift_wt>(r, c) - v2)*second_deriv_scale;
float dxy = (img.at<sift_wt>(r+1 , c+1 ) - img.at<sift_wt>(r+1 , c-1 ) -
img.at<sift_wt>(r-1 , c+1 ) + img.at<sift_wt>(r-1 , c-1 ))*cross_deriv_scale;
float dxs = (next.at<sift_wt>(r, c+1 ) - next.at<sift_wt>(r, c-1 ) -
prev.at<sift_wt>(r, c+1 ) + prev.at<sift_wt>(r, c-1 ))*cross_deriv_scale;
float dys = (next.at<sift_wt>(r+1 , c) - next.at<sift_wt>(r-1 , c) -
prev.at<sift_wt>(r+1 , c) + prev.at<sift_wt>(r-1 , c))*cross_deriv_scale;
Matx33f H(dxx, dxy, dxs,
dxy, dyy, dys,
dxs, dys, dss);
Vec3f X = H.solve(dD, DECOMP_LU);
xi = -X[2 ];
xr = -X[1 ];
xc = -X[0 ];
if ( std ::abs (xi) < 0.5f && std ::abs (xr) < 0.5f && std ::abs (xc) < 0.5f )
break ;
if ( std ::abs (xi) > (float )(INT_MAX/3 ) ||
std ::abs (xr) > (float )(INT_MAX/3 ) ||
std ::abs (xc) > (float )(INT_MAX/3 ) )
return false ;
c += cvRound(xc);
r += cvRound(xr);
layer += cvRound(xi);
if ( layer < 1 || layer > nOctaveLayers ||
c < SIFT_IMG_BORDER || c >= img.cols - SIFT_IMG_BORDER ||
r < SIFT_IMG_BORDER || r >= img.rows - SIFT_IMG_BORDER )
return false ;
}
if ( i >= SIFT_MAX_INTERP_STEPS )
return false ;
关于高斯尺度参数sigma的取值:sigma有两种,一种是相对原始图像的绝对sigma,一种是相对前一层图像的相对sigma 。 step 1: 对于相同octave的图像,设上下层的sigma变化比率为k,设A的绝对sigma为a,A的前一层B的绝对sigma为b,则A相对于B的相对sigma:
c 2 = a 2 − b 2 = ( k 2 − 1 ) b 2
。因此如果取
k = 2 √
,则有c=b. step 2: 对于不同octave的图像,先生成每个octave底层的图像,利用step 1 生成整个octave的图像。设A的绝对sigma为a,A前面的octave B的绝对sigma为b,则有a=2b. 如果
k = 2 1 / i d x
,则octave的图像仅仅需要从上一层octave的第idx张图像C进行直接采样即可,图像C的绝对sigma
c = b ∗ k i d x = 2 b
这样就不需要进行额外的高斯模糊,从而加快速度。
void SIFT::buildGaussianPyramid( const Mat& base, vector <Mat> & pyr, int nOctaves ) const
{
vector <double > sig(nOctaveLayers + 3 );
pyr.resize(nOctaves*(nOctaveLayers + 3 ));
sig[0 ] = sigma;
double k = pow ( 2. , 1. / nOctaveLayers );
for ( int i = 1 ; i < nOctaveLayers + 3 ; i++ )
{
double sig_prev = pow (k, (double )(i-1 ))*sigma;
double sig_total = sig_prev*k;
sig[i] = std ::sqrt (sig_total*sig_total - sig_prev*sig_prev);
}
for ( int o = 0 ; o < nOctaves; o++ )
{
for ( int i = 0 ; i < nOctaveLayers + 3 ; i++ )
{
Mat& dst = pyr[o*(nOctaveLayers + 3 ) + i];
if ( o == 0 && i == 0 )
dst = base;
else if ( i == 0 )
{
const Mat& src = pyr[(o-1 )*(nOctaveLayers + 3 ) + nOctaveLayers];
resize(src, dst, Size(src.cols/2 , src.rows/2 ),
0 , 0 , INTER_NEAREST);
}
else
{
const Mat& src = pyr[o*(nOctaveLayers + 3 ) + i-1 ];
GaussianBlur(src, dst, Size(), sig[i], sig[i]);
}
}
}
}