opencv sift 的实现与原理

reference

opencv sift 原理与原码分析

wikipedia

AI Shack

opencv sift c

一点补充

关于插值：利用了Hessian矩阵中不定点的特性，如果DOG空间某点为极值点，那么该点的Hessian存在特征点（不动点)。可以通过迭代进行插值计算。

//
// Interpolates a scale-space extremum's location and scale to subpixel
// accuracy to form an image feature. Rejects features with low contrast.
// Based on Section 4 of Lowe's paper.
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);
//通过LU分解，求解 HX=dD
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;
}

// ensure convergence of interpolation
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: c2=a2−b2=(k2−1)b2c^2=a^2-b^2=(k^2-1)b^2。因此如果取k=2√k=\sqrt 2，则有c=b.

step 2: 对于不同octave的图像，先生成每个octave底层的图像，利用step 1 生成整个octave的图像。设A的绝对sigma为a,A前面的octave B的绝对sigma为b,则有a=2b. 如果k=21/idxk=2^{1/idx}，则octave的图像仅仅需要从上一层octave的第idx张图像C进行直接采样即可，图像C的绝对sigma c=b∗kidx=2bc=b*k^{idx}=2b 这样就不需要进行额外的高斯模糊，从而加快速度。

//opencv sift source code
void SIFT::buildGaussianPyramid( const Mat& base, vector<Mat>& pyr, int nOctaves ) const
{
vector<double> sig(nOctaveLayers + 3);
pyr.resize(nOctaves*(nOctaveLayers + 3));

// precompute Gaussian sigmas using the following formula:
//  \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2
//sigma为初始值，默认为1.6
sig[0] = sigma;
//nOctaveLayers 默认为3
double k = pow( 2., 1. / nOctaveLayers );
for( int i = 1; i < nOctaveLayers + 3; i++ )
{
double sig_prev = pow(k, (double)(i-1))*sigma;
//绝对sigma
double sig_total = sig_prev*k;
//相对sigma
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;
// base of new octave is halved image from end of previous octave
else if( i == 0 )
{
const Mat& src = pyr[(o-1)*(nOctaveLayers + 3) + nOctaveLayers];
//src 是上一个octave的第nOctaveLayers张图像，恰好满足2倍sigma的要求。这样只需要resize就可以得到下一个octave的图像了。
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]);
}
}
}
}