SIFT特征（三）

上一章节中，我们探讨了如何检测尺度空间上的极值点，以及如何排除掉不稳定的极值点。为了使得SIFT具有旋转不变性，我们需要利用关键点的邻域信息为每一个关键点分配一个方向，使得其具有局部结构的稳定方向。

5.关键点方向赋值
5.1 梯度方向及幅值的计算
对于DoG金字塔中检测出的稳定关键点，取关键点邻域$3*\sigma$ 大小的窗口内的像素点，并分别计算每个像素点的梯度方向和大小，公式如下，

$$m(x,y)=\sqrt{(L(x+1,y)-L(x-1,y))^2+(L(x,y+1)-L(x,y-1))^2}$$
$$\theta(x,y)=tan^{-1}((L(x,y+1)-L(x,y-1))/(L(x+1,y)-L(x-1,y)))$$
$L$ 为DoG空间，这里邻域高斯函数的参数为$\sigma=1.5*\sigma_{oct}$， $\sigma_{oct}$ 为关键点的尺度大小。按照采样 $3*\sigma$ 原则，取领域半径为 $3*\sigma=3*1.5*\sigma_{oct}$。见下图左侧图所示，蓝色圈表示高斯作用区域，黑色箭头指示每个对应像素点的梯度。

5.2 计算梯度直方图
上一步中计算了邻域内各个像素点的梯度大小和方向，接下来使用直方图对这些邻域像素点的方向和大小进行统计。直方图按照0-360度的方向范围划分为36个柱（bins），即每个柱的度数范围为10度。依次遍历邻域内每个像素点，按照像素点梯度方向找到该像素点所对应的柱的编号，并将像素点的梯度大小与权重的乘积累加到所对应的直方图柱中，这里的权重是指采用了参数为$\sigma$高斯函数对邻域像素梯度值进行处理，这样加权处理后，距离关键点较远的边缘像素点的贡献值将会降低，而距离关键点距离较近的像素点则具有较多的贡献值。
上图右侧图示意将关键点邻域像素点根据其梯度方向和大小统计后的直方图。我们可以看出，在统计直方图中，会存在一个峰值方向，那么我们将这个峰值方向代表该关键点的主方向。为了增强匹配的鲁棒性，另外保留了直方图柱大于主方向峰值大小80%的方向，作为关键点的辅方向。在这种情况下，在相同的位置和尺度将会有多个关键点被创建但方向不同，其方向分别表示主方向与各个辅方向。如下图表示关键点检测示意图。

---

5.3 C++源码解读

/*计算关键点邻域的梯度直方图
* pt:关键点xy坐标
*radius: 3*1.5* scale，邻域窗口半径
* sigma：1.5*scale，高斯函数sigma值
*hist: 梯度方向直方图
* n: 直方图bin数量
*/
static float calcOrientationHist( const Mat& img, Point pt, int radius, float sigma, float* hist, int n )
{
    int i, j, k, len = (radius*2+1)*(radius*2+1);

    float expf_scale = -1.f/(2.f * sigma * sigma);
    AutoBuffer<float> buf(len*4 + n+4);
    float *X = buf, *Y = X + len, *Mag = X, *Ori = Y + len, *W = Ori + len;
    float* temphist = W + len + 2;

    for( i = 0; i < n; i++ )
        temphist[i] = 0.f;

    //遍历邻域内所有像素点统计直方图
    for( i = -radius, k = 0; i <= radius; i++ )
    {
        int y = pt.y + i;
        if( y <= 0 || y >= img.rows - 1 )
            continue;
        for( j = -radius; j <= radius; j++ )
        {
            int x = pt.x + j;
            if( x <= 0 || x >= img.cols - 1 )
                continue;

            //计算当前像素点xy方向的梯度
            float dx = (float)(img.at<sift_wt>(y, x+1) - img.at<sift_wt>(y, x-1));
            float dy = (float)(img.at<sift_wt>(y-1, x) - img.at<sift_wt>(y+1, x));

            //X存储x方向梯度大小；Y存储y方向梯度大小；W存放像素点权重
            X[k] = dx; Y[k] = dy; W[k] = (i*i + j*j)*expf_scale;
            k++;
        }
    }

    //直方图统计总的像素点个数
    len = k;

    // 计算邻域内像素点的梯度大小（magnitude）、方向（fastAtan2）和权值（exp）
    //计算矩阵元素W的自然e的W次幂。i.e. W=exp(W)。值为：exp((i^2+j^2)/2*sigma^2)，距离关键点越远，权重越小
    exp(W, W, len);
    //计算梯度方向，存放在Ori中
    fastAtan2(Y, X, Ori, len, true);
    //计算梯度幅值，存放在Mag中。i.e. mag（I） = \sqrt(X^2+Y^2), X和Y 是xy方向的梯度大小，可以根据sobel等算法计算而来
    magnitude(X, Y, Mag, len);

    //计算每个像素点所在bin的ID，并将值梯度大小*权重放在所对应的bin中
    for( k = 0; k < len; k++ )
    {
        int bin = cvRound((n/360.f)*Ori[k]);
        if( bin >= n )
            bin -= n;
        if( bin < 0 )
            bin += n;
        temphist[bin] += W[k]*Mag[k];
    }

    // 平滑直方图（取距离为3的直方图柱插值）
    temphist[-1] = temphist[n-1];
    temphist[-2] = temphist[n-2];
    temphist[n] = temphist[0];
    temphist[n+1] = temphist[1];
    for( i = 0; i < n; i++ )
    {
        hist[i] = (temphist[i-2] + temphist[i+2])*(1.f/16.f) +
            (temphist[i-1] + temphist[i+1])*(4.f/16.f) +
            temphist[i]*(6.f/16.f);
    }

    //计算关键点的主方向
    float maxval = hist[0];
    for( i = 1; i < n; i++ )
        maxval = std::max(maxval, hist[i]);

    return maxval;
}