sift特征--关键点搜索与定位

     前面介绍了如何生成高斯图像金字塔，并计算了每组图像的高斯差分图像。现在介绍如何进行关键点搜索与定位（都在灰度图上搞的）。

一、极值点计算
   关键点是由DOG空间的局部极值点组成的。为了寻找DoG函数的极值点，每一个像素点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。



每一幅高斯差分图像中的一个像素点，要和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。
具体算法：
(1)   遍历不同分辨率和尺度高斯差分金字塔，获取特定分辨率和尺度下的高斯差分图像Io,s
(2)   同组同层比较，对Io,s内的每个像素Pcenter，首先比较Pcenter与阈值PRE_COLOR_THRES，若不满足条件，则返回false，即该点不是极值点，转到(4)；否则，将Pcenter与其周围8邻域像素点比较，只要该邻域内存在相对于Pcenter的最大和最小值像素点，则认为Pcenter是该层的极值点，返回true，转到(3)
(3)   同组不同层比较，在Io,s上下两层高斯差分图像Io,s-1和Io,s+1中，分别将Pcenter同其对应的上下层邻域进行比较，如果在上下两层邻域内均存在相对于Pcenter的最大最小值，则认为Pcenter是相邻层的极值点，返回true，否则就不是极值点，返回false，转到(4)
(4)   对于返回值为true的，将Pcenter加入极值点坐标集合ret中
 
代码如下，清华小哥用到了c++的新特性，lamda函数，略高端：vector<Coor> ExtremaDetector::get_raw_extrema() const {
vector<Coor> ret;
int npyramid = dog.noctave, nscale = dog.nscale;
REP(i, npyramid)
REPL(j, 1, nscale - 2) {
auto& now = dog.dogs[i][j]; //dogs有n组，每组有s-1层尺度差分图像
int w = now.width(), h = now.height();

auto v = get_local_raw_extrema(i, j); //返回尺度差分图像邻域26个点内满足参考阈值的中心检测点坐标集合
for (auto& c : v) {
ret.emplace_back((float)c.x / w * dog.origw, //w,h表示当前尺度图像的宽和高，orignw,orignh表示原始图像(第1组第1层)的宽和高
(float)c.y / h * dog.origh); //将检测到的特征点坐标统一到原始图像上
}
}
return ret;
}

vector<Coor> ExtremaDetector::get_local_raw_extrema(
int pyr_id, int scale_id) const {
vector<Coor> ret;

const Mat32f& now(dog.dogs[pyr_id][scale_id]); //获取特定分辨率和尺度下的高斯差分图像
int w = now.width(), h = now.height();

//lamda函数
auto is_extrema = [this, &now, pyr_id, scale_id](int r, int c) { //is_extrema为函数名，中括号为了共用父函数的变量，圆括号为is_extrema的入参
float center = now.at(r, c);
if (center < PRE_COLOR_THRES) // initial color is less than thres
return false;

bool max = true, min = true;
float cmp1 = center - JUDGE_EXTREMA_DIFF_THRES,
cmp2 = center + JUDGE_EXTREMA_DIFF_THRES;
// try same scale
REPL(di, -1, 2) REPL(dj, -1, 2) { //同组同层(同一尺度) 比较，8邻域，(0,0)被忽略
if (!di && !dj) continue;
float newval = now.at(r + di, c + dj);
if (newval >= cmp1) max = false; //newval对于cmp1和cmp2不是同一个值，只要是该循环内都行，我快疯了$_$
if (newval <= cmp2) min = false; //对于该中心检测点，只要在它的邻域内存在满足阈值的最大最小点，就记录该中心监测点，下同
if (!max && !min) return false;
}

// try adjencent scale
for (int ds = -1; ds < 2; ds += 2) { //同组不同层比较，ds表示相对层数，nl即不同尺度图像的层数索引
int nl = scale_id + ds;
auto& mat = this->dog.dogs[pyr_id][nl];

REPL(di, -1, 2) {
const float* p = mat.ptr(r + di) + c - 1; //p为中心点监测点在相邻层的对应点
REP(i, 3) {
float newval = p[i];
if (newval >= cmp1) max = false;
if (newval <= cmp2) min = false;
if (!max && !min) return false;
}
}
}
return true;
};

REPL(i, 1, h - 1) REPL(j, 1, w - 1)
if (is_extrema(i, j))
ret.emplace_back(j, i);
return ret;
}

   有一个问题是到底要在多少个尺度中寻找极值点， 即如何确定 s 值(组内层数)。实验表明， s 取 3 是较好的选择。如果 s = 3，则需要 5 幅高斯差分图像才可以。这里的计算是高效的，因为大多数情况下，只需要几步比较，就可以排除一个像素点，认为它不是极值。
另外，在极值比较的过程中，每一组图像的首末两层是无法进行极值比较的，为了满足尺度变化的连续性，我们在每一组图像的顶层继续用高斯模糊生成了3幅图像，高斯金字塔有每组S+3层图像。DOG金字塔每组有S+2层图像。也就是说，我们只牺牲了-1组的第0层和第N组的最高层。
二、抽取稳定的关键点
   由于DoG值对噪声和边缘较敏感,因此,在上面DoG尺度空间中检测到局部极值点还要经过进一步的检验才能精确定位为特征点。上一步已经求出了极值点，现在要对这些极值点进行筛选，去除不稳定的点，以增强特征点匹配时的稳定性、提高抗噪声能力。不稳定的点包括低对比度的点和边缘上的点。同时，由于在金子塔中存在降采样的图像，在这些图像中提取的极值点在原始输入图像中到底在什么位置，也是一个问题。下面将提出上面两个问题的解决方案。
 
1.      去除低对比度点
   低对比度的极值点，是那些对噪声敏感的候选点 ，需要剔除。遍历每个候选极值点sp，利用插值法重新计算坐标，并重复多次后，再将新坐标的像素与阈值比较，决定去留。对于保留的极值点，更新其坐标，进入边缘响应剔除步骤。
   插值的计算采用对尺度空间DoG函数进行曲线拟合，DoG函数在尺度空间的Taylor展开式，即在某极值点A对D(x,y,σ)进行泰勒展开：


         (1)    这里X=(x,y, σ)是最终插值点到点A的偏移量，对上式求X的偏导数，并令之为零，求X


                                                                                (2)
     Sift原著作者认为如果，偏移量中任何一个分量大于0.5，则这个极值点和另一个采样点（图像中的另一个像素）离得更近，需要采用插值法求得极值点位置的估计值(即分量大于0.5，offset不可忽略~~)，同时也可以将x待回原式求得D(x)与阈值0.03比较去除低对比度点：


                                                            (3)
     上面的式子中，二阶导数用Hessian矩阵，一阶导数用梯度，它们均通过相邻像素的差值计算。这里的Hessian矩阵是由三元函数的二阶偏导数构成的方阵，描述了DOG函数的局部曲率：




                                                                    (4)

具体算法流程：
    遍历每个极值点，利用(4)式计算新的插值坐标偏移量iter_offset。若满足OFFSET_THRES阈值要求，则更新极值坐标，并将新坐标代入(3)式，求得的结果dextr与阈值CONTRAST_THRES比较，若满足要求，则更新极值点的真实坐标sp->real_coor
vector<SSPoint> ExtremaDetector::get_extrema() const {
TotalTimer tm("extrema");
int npyramid = dog.noctave, nscale = dog.nscale;
vector<SSPoint> ret;
#pragma omp parallel for schedule(dynamic) //不同分辨率并行处理
REP(i, npyramid)
REPL(j, 1, nscale - 2) {
auto v = get_local_raw_extrema(i, j); //获取极大值
//print_debug("raw extrema count: %lu\n", v.size());
for (auto& c : v) {
SSPoint sp;
sp.coor = c;
sp.pyr_id = i;
sp.scale_id = j;
bool succ = calc_kp_offset(&sp); //去除低对比度
if (not succ) continue;
auto& img = dog.dogs[i][sp.scale_id];
succ = not is_edge_response(sp.coor, img); //去除边缘响应
if (not succ) continue;

#pragma omp critical
ret.emplace_back(sp);
}
}
return ret;
}

/*去除对比度较低的不稳定极值点*/
bool ExtremaDetector::calc_kp_offset(SSPoint* sp) const {
TotalTimer tm("offset");
auto& now_pyramid = dog.dogs[sp->pyr_id]; //now_pyramid同一分辨率图数组，包含一组尺度图像
auto& now_img = now_pyramid[sp->scale_id]; //now_img当前尺度图像(差分+模糊)
int w = now_img.width(), h = now_img.height();
int nscale = dog.nscale;

Vec offset, delta; // partial(d) / partial(offset)
int nowx = sp->coor.x, nowy = sp->coor.y, nows = sp->scale_id; //sp表示局部极大值的候选点
int niter = 0;
for(;niter < CALC_OFFSET_DEPTH; ++niter) { //对某个特征点重复进行的修正次数，而非遍历什么图像...
if (!between(nowx, 1, w - 1) ||
!between(nowy, 1, h - 1) ||
!between(nows, 1, nscale - 2))
return false;

auto iter_offset = calc_kp_offset_iter(
now_pyramid, nowx, nowy, nows);
offset = iter_offset.first,
delta = iter_offset.second;
if (offset.get_abs_max() < OFFSET_THRES) // found，If the offset X(x,y,z) is larger than 0.5 in any dimension, then it means that the extremum lies closer to a different sample point，由于是任何尺度，故直接求max了
break;

nowx += round(offset.x); //calc_kp_offset_iter得到的坐标均大于OFFSET_THRES,需要进行插值估计，重新计算
nowy += round(offset.y);
nows += round(offset.z);
}
if (niter == CALC_OFFSET_DEPTH) return false;

double dextr = offset.dot(delta); // calc D(x~)，delta*X
dextr = now_pyramid[nows].at(nowy, nowx) + dextr / 2;
// contrast too low
if (dextr < CONTRAST_THRES) return false;

// update the point
sp->coor = Coor(nowx, nowy);
sp->scale_id = nows;
sp->scale_factor = GAUSS_SIGMA * pow(
SCALE_FACTOR, ((double)nows + offset.z) / nscale); //该delta为尺度图像进行高斯滤波的窗口大小时的delta
// accurate real-value coor
sp->real_coor = Vec2D(
((double)nowx + offset.x) / w * dog.origw,
((double)nowy + offset.y) / h * dog.origh);
return true;
}

/*返回精细化后局部极值的具体坐标及delta，供calc_kp_offset进一步阈值挑选*/
std::pair<Vec, Vec> ExtremaDetector::calc_kp_offset_iter(
const DOGSpace::DOG& now_pyramid,
int x , int y, int s) const {
Vec offset = Vec::get_zero();
Vec delta;
double dxx, dyy, dss, dxy, dys, dsx;

auto& now_scale = now_pyramid[s];
#define D(x, y, s) now_pyramid[s].at(y, x) //D和DS的区别在于，D可以随意选择尺度索引S
#define DS(x, y) now_scale.at(y, x)
float val = DS(x, y);

//计算梯度---> D关于X的一阶偏导数(负的)
delta.x = (DS(x + 1, y) - DS(x - 1, y)) / 2;
delta.y = (DS(x, y + 1) - DS(x, y - 1)) / 2;
delta.z = (D(x, y, s + 1) - D(x, y, s - 1)) / 2; //delta.z是指前一层与后一层两层差分图像的之间的差值

//计算Hessian矩阵，描述函数的局部曲率---->表示D关于X的二阶偏导数
dxx = DS(x + 1, y) + DS(x - 1, y) - val - val;
dyy = DS(x, y + 1) + DS(x, y - 1) - val - val;
dss = D(x, y, s + 1) + D(x, y, s - 1) - val - val;

dxy = (DS(x + 1, y + 1) - DS(x + 1, y - 1) - DS(x - 1, y + 1) + DS(x - 1, y - 1)) / 4;
dys = (D(x, y + 1, s + 1) - D(x, y - 1, s + 1) - D(x, y + 1, s - 1) + D(x, y - 1, s - 1)) / 4;
dsx = (D(x + 1, y, s + 1) - D(x - 1, y, s + 1) - D(x + 1, y, s - 1) + D(x - 1, y, s - 1)) / 4;
#undef D
#undef DS

Matrix m(3, 3); //m表示D关于X的二阶偏导数
m.at(0, 0) = dxx; m.at(1, 1) = dyy; m.at(2, 2) = dss;
m.at(0, 1) = m.at(1, 0) = dxy;
m.at(0, 2) = m.at(2, 0) = dsx;
m.at(1, 2) = m.at(2, 1) = dys;

Matrix pdpx(3, 1); // delta = dD / dx
delta.write_to(pdpx.ptr());

Matrix inv;
if (not m.inverse(inv)) { // pseudo inverse is slow，对m求逆矩阵，用SVD方法pseudo_inverse
inv = m.pseudo_inverse();
}
auto prod = inv.prod(pdpx); //矩阵相乘，x,y二维坐标，z表示尺度坐标
offset = Vec(prod.ptr());
return {offset, delta};
}

2.      去除边缘响应点
    由于DoG函数在图像边缘有较强的边缘响应，而边缘上的极值点抗噪性较差，因此我们还需要排除边缘响应。
    我们知道曲面上每个点（非平点）都有两个主方向，并且沿这两个主方向的法曲率（即两个主曲率）分别是曲面在该点法曲率的最大值和最小值。对于边缘上的点，沿垂直于边缘的方向上，法曲率最大，而沿边缘的方向上，法曲率最小。因此对于分布在边缘上附近的极值点，它们的法曲率最大值和最小值之比（即两个主曲率之比），一般情况下要比非边缘点的比值大。根据这种思想，我们可以设一个比值的阈值，当比值大于这个阈值就认为极值点在边缘上。
    DOG函数主曲率可以通过计算在该点位置尺度的2×2的Hessian矩阵得到，导数由采样点相邻差来估计：


    这里Dxx表示DOG金字塔中某一尺度的图像x方向求导两次，微分可以通过计算邻近点的差值来近似计算。
    又因为DOG的主曲率和H的特征值成正比，我们只需要计算 H 的较大特征值与较小特征值的比例即可。设α 是较大的特征值， β 是较小的特征值，由矩阵性质知：



    其中用到了矩阵的迹和行列式。通常这里的行列式不会是负值，如果出现负值的情况，即两个主曲率不同号，我们将丢弃这个点，不将其视为极值点。设 r=α/ β ，我们可得：



      这里当 r≥1，(r+1)2/r是r的单调递增函数，因此要计算主曲率的比值（即 r）是否在某阈值之下，为了避免求H的特征值，只需要判断上式左边的项是否在阈值之下即可，通常r取10。
 
具体步骤：
  对于每个候选特征点，计算其Hessian矩阵，然后计算该矩阵的迹和行列式的比值，如果该比值大于阈值则保留该候选关键点，否则剔除。
/*去除边缘响应*/
bool ExtremaDetector::is_edge_response(Coor coor, const Mat32f& img) const {
float dxx, dxy, dyy;
int x = coor.x, y = coor.y;
float val = img.at(y, x);

//计算Hessian矩阵
dxx = img.at(y, x + 1) + img.at(y, x - 1) - val - val;
dyy = img.at(y + 1, x) + img.at(y - 1, x) - val - val;
dxy = (img.at(y + 1, x + 1) + img.at(y - 1, x - 1) -
img.at(y + 1, x - 1) - img.at(y - 1, x + 1)) / 4;
float det = dxx * dyy - dxy * dxy;
if (det <= 0) return true;
float tr2 = sqr(dxx + dyy);

// Calculate principal curvature by hessian
/*极值点分布在边缘上，其主曲率比值要比非边缘点的比值大，因此只需计算主曲率比值即可，
* 而由于主曲率与海森矩阵的特征值成正比，咱只需计算海森矩阵的特征值之比即可*/
if (tr2 / det < sqr(EDGE_RATIO + 1) / EDGE_RATIO) return false;
return true;
}

极大值点结果（未去除低对比度和边缘响应）


精确定位特征点处理效果：



    注意这里的特征点来自于不同分辨率的不同尺度的高斯差分图像，下图展示的是经过统一坐标大小后的结果。