SIFT特征（二）

在上一章节中，我们介绍了如何对图像模拟不同尺度空间，以在尺度空间下检测稳定特征点。构造尺度空间的方法是通过构造高斯金字塔和DoG实现的。那么在本章节中，我们将讨论如何在构造好的DoG空间内检测关键点，以及筛选出不稳定的关键点。

4.检测并筛选关键点（keypoints）
4.1 检测尺度空间极值点
计算DoG空间的极值点（局部最大值或最小），是通过这样的方法实现的：每一个像素点需要与同一层内相邻的8个像素点以及相邻上下两层的$9*2=18$ 个，共26个像素点比较大小。如果当前像素点比这26个像素点都要大或者都要小，则当前像素点为极值点。如下图所示，其中x号像素表示当前像素，绿色圈表示像素为相邻比较像素，

4.2确定检测点的精确位置
实际上，经过这样步骤检测的极值点是离散的，可能并不是真实的极值点，这可以从下面这个示意图看出，真正的极值点是位于检测的离散极值点附近的。所以为了得到关键点的准确位置，Lowe是通过对关键点附近的离散点进行拟合，然后求极值，从而得到关键点的准确位置的。

为了说明对离散极值点拟合逼近真实极值点的过程，我们首先举出一个一维函数插值的例子。
假设有，$f(-1)=1,f(0)=6,f(1)=5$， 求$f(x)$ 在区间 [-1,1]上的最大值。
$f(x)$ 在$f(0)$ 附近的泰勒展开式为，

$$f(x)\approx f(0)+f'(0)x+\frac{f''(0)}{2}x^2$$
另外有， $f(x)$ 在 $x$ 处的导数可以写成离散的形式为，$f'(x)=\frac{f(x+1)-f(x-1)}{2}$ ，二阶导数写成离散形式为， $f''(x)=f(x+1)+f(x-1)-2f(x)$
所以有，
$$f(x)\approx 6-2x-3x^2$$
对函数 $f(x)$ 求导得到极值点的坐标为 $\hat{x}=\frac{1}{3}$， 则极值点的坐标为 $f(\hat{x})=6\frac{1}{3}$

同样的道理，我们对检测出的离散关键点进行拟合求极值，假设在尺度$\sigma$ 的图像空间$D(x,y)$ 处检测到一个极值点（关键点），那么这个极值点的空间坐标为$(x,y,\sigma)$ ，通过上面分析我们知道，真实的极值点有可能落在这个离散极值点的附近。设$(x,y,\sigma)$ 附近的点相对于该极值点的偏移量为$(\Delta x,\Delta y, \Delta \sigma)$,那么$D(\Delta x,\Delta y, \Delta \sigma)$ 可以表示为点$(x,y,\sigma)$ 处的泰勒展开式，

$$D(X)=D+\frac{\partial D^T}{\partial X}X+\frac{1}{2}X^{T}\frac{\partial ^2D}{\partial X^2}X$$
$X$ 表示偏移量$(\Delta x,\Delta y, \Delta \sigma)$，对上面的式子求一阶导数并等于0，则得到
$$X=-\frac{\partial^2D^{-1}}{\partial x^2}{\frac{\partial D}{\partial x}}$$
多次迭代，每次更新极值点的位置， Lowe算法中最多迭代5次，最终得到候选极值点的精确位置 $\hat{x}$。
$$x\leftarrow x+X$$
将极值点的精确位置 $\hat{x}$ 带入 $D$，得到极值点的像素值。
如果$|D(\hat{x})|$ 的值过小（Lowe算法中，设置阈值为0.03），说明该点处像素值很容易受到噪声的干扰而变得不稳定，所以需要将该点剔除掉。

4.3 消除边缘响应
DoG对图像中的边缘具有较强的响应，而一旦特征点落在边缘上，这些点就是不稳定的检测点，所以需要将这些边缘点去掉。
而在横跨边缘处，DoG往往具有较大的主曲率，而在垂直边缘的方向具有较小的主曲率，所以我们可以考虑，如果检查最大主曲率与最小主曲率的比值，如果比值大于某一阈值，就认为是边缘点，并剔除该关键点，就可以将边缘不稳定的关键点排除掉。
海塞矩阵的特征值是和D的主曲率是成正比的，即计算主曲率可以转换为计算海塞矩阵。

$$H(x,y)= \begin{bmatrix} D_xx & D_xy \\ D_xy & D_yy \\ \end{bmatrix}$$
$D$ 可以通过邻近像素点求差分得到。设H的最大特征值和最小特征值分别为$\alpha$和 $\beta$，则还塞矩阵的迹和行列式表示为，
$$Tr(H) = D_{xx}+D_{yy}=\alpha+\beta$$
$$Det(H)=D_{xx}D_{yy}-(D_{xy})^2=\alpha\beta$$
如果 $\gamma$ 表示最大特征值和最小特征值之间的比例，即 $\alpha=\gamma\beta$，那么有，
$$\frac{Tr(H)^2}{Det(H)}=\frac{(\alpha+\beta)^2}{\alpha\beta}=\frac{(\gamma+1)^2}{\gamma}$$
当两个特征值相等时， $\frac{(\gamma+1)^2}{\gamma}$ 有最小值，而随着两者比例 $\gamma$ 的增加， $\frac{(\gamma+1)^2}{\gamma}$ 也增加。所以想要确定D的主曲率是否大于某一阈值 $\hat{\gamma}$， 只需要检查下面的式子是否成立，
$$\frac{Tr(H)^2}{Det(H)}>\frac{(\hat{\gamma}+1)^2}{\hat{\gamma}}$$
在Lowe算法中，给定的阈值 $\hat{\gamma}$是10，也就是说当最大主曲率与最小主曲率之间的比值大于10时，该关键点被删除。

经过上面三个步骤，即首先通过比较像素点与同层8邻域像素点和相邻上下两层像素点确定候选关键点；然后通过曲线拟合确定关键点的准确位置，并删除对比度较低的极值点；最后通过判断最大主曲率与最小主曲率之间的比值，删除边缘DoG相应较强的关键点。最终得到的关键点用于SIFT特征计算。

---

4.4 C++源码解读
1. 曲线拟合精确确定关键点的坐标

    //迭代插值过程确定关键点精确坐标
    for( ; i < SIFT_MAX_INTERP_STEPS; i++ )
    {
        //octv当前关键点的组数；layer当前关键点的层数；nOctaveLayers每组的总层数
        //dog_pyr:DoG金字塔
        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];

        //如果求得的偏移极值点大于0.5，则认为真实的极值点是位于其他的采样点的，此时需要更新极值点的坐标
        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;

2.判断最大主曲率与最小主曲率比值来判断是否是边缘

{       
        //octv当前关键点的组数；layer当前关键点的层数；nOctaveLayers每组的总层数
        //dog_pyr:DoG金字塔
        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];

        Matx31f 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 t = dD.dot(Matx31f(xc, xr, xi));

        contr = img.at<sift_wt>(r, c)*img_scale + t * 0.5f;
        if( std::abs( contr ) * nOctaveLayers < contrastThreshold )
            return false;

        // 利用海塞矩阵的迹和行列式计算D主曲率的比值
        float v2 = img.at<sift_wt>(r, c)*2.f;
        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 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 tr = dxx + dyy;//迹
        float det = dxx * dyy - dxy * dxy;//行列式

        //判断主曲率的比例是否小于阈值edgeThreshold，如果不是，则认为是边缘关键点，删除该关键点
        if( det <= 0 || tr*tr*edgeThreshold >= (edgeThreshold + 1)*(edgeThreshold + 1)*det )
            return false;
    }